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Abstract 

Random  phase  screens  are  essential  elements  of  simulating  light  propagation 
through  turbulent  media.  In  order  to  be  effective,  they  must  accurately  reflect  theory 
and  be  implementable  by  the  user.  This  document  explains  and  evaluates  three 
methods  of  generating  random  phase  screens:  using  a  Fourier  series  upon  a  polar 
frequency  grid  with  logarithmic  spacing;  using  the  fast  Fourier  transform,  with  its 
cartesian  frequency  grid;  and  using  Zernike  polynomials.  It  provides  a  comparison  of 
the  polar  Fourier  series  technique  with  the  two  more  common  techniques  (fast  Fourier 
transform  and  Zernike),  with  the  end  result  of  giving  the  users  enough  information  to 
choose  which  method  best  fits  their  needs.  The  evaluation  criteria  used  are  generation 
time  (usability)  and  phase  structure  function  (accuracy). 
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Polar  Phase  Screens: 


A  Comparison  with  Other  Methods 

of 

Random  Phase  Screen  Generation 

I.  Overview 

Random  phase  screens  are  key  elements  of  simulating  optical  atmospheric  tur¬ 
bulence.  There  several  methods  of  generating  random  phase  screens,  three 
of  which  (Zernike,  Fast  Fourier  Transform,  and  Polar  Fourier  Series)  are  compared 
herein.  This  is  the  first  organized  analysis  of  the  benefits  and  drawbacks  of  each 
method  using  the  structure  function  as  a  measure  of  accuracy  to  theory.  The  phase 
screen  generation  time  is  used  as  a  measure  of  usability.  Together,  the  structure  func¬ 
tion  and  the  generation  time  allow  one  to  choose  the  most  effective  method  for  a  given 
application. 

1.1  Why  Random  Phase  Screens ? 

Light  travelling  through  the  atmosphere  (or  any  turbulent  medium)  is  affected 
by  the  random  fluctuations  of  the  index  of  refraction  (IOR).  Coherent  masses  of  air, 
called  eddies,  are  characterized  by  indices  of  refraction  determined  by  their  partic¬ 
ular  temperature,  pressure  and  humidity.  The  effect  of  these  eddies  upon  the  light 
incident  upon  them  is  observable  as  distortion,  or  even  loss,  of  the  signal  incident 
upon  the  detector.  If  the  optical  path  were  non-random,  the  system  could  be  mod¬ 
elled  and  each  distortion  compensated  for,  allowing  the  perfect  reconstruction  of  the 
source.  However,  as  the  optical  path  is  random,  complex  algorithms  are  required  to 
statistically  model  the  random  data. 

Random  phase  screens  are  used  to  model  the  random  IOR  of  a  turbulent 
medium.  They  do  this  by  adding  a  correlated  random  phase  to  the  propagating 
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field,  resulting  in  a  ‘detected’  image  similar  to  that  which  would  be  the  result  of  ac¬ 
tual  propagation.  More  than  one  screen  can  be  used  to  simulate  a  complete  optical 
path;  and,  because  each  phase  screen  compresses  the  effects  of  propagation  through  a 
three-dimensional  volume  into  a  two-dimensional  ‘screen’,  they  can  be  placed  at  any 
location  along  the  propagation  path. 

There  are  several  reasons  for  using  more  than  one  phase  screen  to  simulate  a 
propagation.  For  instance,  consider  a  telescope  viewing  an  object  in  space.  Using 
only  a  single  phase  screen  and  placing  it  in  the  lens  of  the  receiving  aperture  would 
entirely  negate  the  effect  of  scintillation.  Scintillation,  or  fluctuation  of  signal  power,  is 
effected  by  a  phase  screens  located  away  from  the  aperture.  Also,  typical  Fourier  phase 
screen  generation  techniques  have  difficulty  representing  low  frequency  abberations 
accurately.  However,  placing  the  screen  near  the  source  or  even  equidistant  between 
the  source  and  detector,  fails  to  capture  the  fact  that  most  of  the  turbulence  is  in 
the  first  few  100  meters  of  above  an  earth  based  telescope.  These  difficulties  can  be 
overcome  by  using  multiple  screens  placed  at  strategic  locations  along  the  propagation 
path.  Because  the  low  frequency  tip  and  tilt  abberations  account  for  approximately 
80%  of  the  system’s  aberrant  power,  artificial  limitations  upon  them  can  greatly 
decrease  the  efficacy  of  the  model. 

The  use  of  multiple  phase  screens  has  some  drawbacks.  Each  screen  takes  a 
certain  amount  of  time  to  generate,  depending  upon  its  size  and  the  technique  used. 
Some  applications  call  for  the  generation  of  extremely  large  phase  screens  (e.g.,the 
implementation  of  Taylor’s  frozen  turbulence  hypothesis  to  simulate  time  dependent 
fluctuations  [1,  59]).  Using  current  techniques,  the  generation  of  multiple  large  phase 
screens  can  take  a  prohibitive  amount  of  time. 

As  phase  screens  are  only  a  tool  used  to  test  the  efficacy  of  what  are  often  cum¬ 
bersome  algorithms,  it  is  important  that  they  be  as  efficient  as  possible.  Minimizing 
the  generation  time  and  accuracy  of  phase  screens  maximizes  their  usefulness  to  the 
researcher. 
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1.2  Thesis  Overview 


Chapter  two  describes  the  theoretical  and  mathematical  basis  for  phase  screens 
and  their  analysis. 

Chapter  three  describes  particular  methods  used  to  generate  and  evaluate  phase 
screens,  including  MATLAB®  commands  and  mathematical  formulae. 

Chapter  four  gives  the  results  of  the  comparative  analysis  performed  using  the 
method  of  generation  and  evaluation  outlined  in  chapter  three. 

Chapter  five  presents  the  conclusions  expands  to  possible  applications  and  rec¬ 
ommendations. 

1 . 3  Contributions 

This  document  gives  an  technical  overview  and  evaluation  of  three  methods  of 
generating  random  phase  screens:  Zernike,  Fast  Fourier  Transform,  and  Polar  Fourier 
Series,  ft  provides  information  for  someone  using  random  phase  screens  to  decide 
which  method  best  fits  their  needs. 
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II.  Atmospheric,  Mathematical  and  Simulation  Theory 


Phase  screens  are  valuable  tools  in  simulating  turbulent  media  (such  as  the  earth’s 
atmosphere)  and  their  effects  upon  the  propagation  of  light.  This  chapter  will 
highlight  their  significance  after  presenting  necessary  background  material.  Section 
2.1  gives  a  rough  overview  of  atmospheric  turbulence  and  how  it  affects  the  propa¬ 
gation  of  light.  Section  2.2  discusses  random  variables  and  some  of  their  statistics. 
Random  variable  are  used  extensively  in  optical  atmospheric  modelling  and  2.2  gives 
and  explanation  of  some  of  the  basic  tools  used.  Section  2.3  expands  upon  those  basic 
tools  by  development  of  the  structure  function.  Section  2.4  covers  the  most  common 
models  of  the  atmospheric  power  spectrum,  including  Kolmogorov’s.  Section  2.5  uses 
the  structure  function  and  power  spectrum  to  define  key  atmospheric  descriptors,  in¬ 
cluding  the  Fried  radius.  Section  2.6  gives  an  overview  of  the  types  of  phase  screens 
focused  upon  in  this  thesis.  Section  2.7  defines  Zernike  polynomials,  their  covariance 
and  useful  related  details.  Section  2.8  covers  the  Fourier  transform  as  it  is  used  in 
creating  phase  screens.  It  also  covers  the  distinctions  between  the  three  methods. 

2.1  Atmospheric  Turbulence 

Kolmogorov  modelled  the  atmosphere  as  a  viscous  fluid,  subject  to  conditions 
of  turbulent  and  laminar  flow  [1,  Ch.  3].  The  nature  of  the  flow  can  be  characterized 
by  a  turbulent  intensity  using  the  Reynolds  Number  [9,  p.  58],  which  increases  with 
the  degree  of  turbulence.  Turbulence  in  the  atmosphere  is  caused  by  local  unstable  air 
masses  resulting  from  convection  and  wind  shear  in  the  atmosphere.  These  air  masses, 
called  eddies  range  in  size  from  small  (on  the  order  of  millimeters  or  centimeters)  to 
large  (on  the  order  of  100  meters).  The  smallest  an  eddy  can  be  and  still  retain  its 
own  distinct  refraction  characteristics  is  called  the  inner  scale,  denoted  Iq.  Thus  air 
masses  smaller  than  l0  are  seen  only  as  elements  of  larger  eddies.  The  largest  and 
eddy  can  be  and  still  maintain  viscosity,  seen  as  a  correlated  index  of  refraction  (IOR), 
is  called  the  outer  scale,  denoted  L0.  The  range  of  eddy  sizes  Iq  to  L0  is  called  the 
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inertial  subrange,  within  which  there  is  a  continuous  transfer  of  energy  as  large  eddies 
break  down  into  smaller  ones  and  small  eddies  eventually  dissipate  as  heat. 

The  defining  characteristic  of  an  eddy  is  that,  within  itself,  each  eddy  is  homo¬ 
geneous  and  isotropic.  The  inertial  subrange  is  the  complete  range  of  eddy  sizes  of 
interest.  Air  masses  smaller  than  Iq  are  assumed  not  to  exist  and  air  masses  larger 
than  L0  are  no  longer  assumed  homogeneous  or  isotropic.  It  should  be  noted  that  the 
inner  scale  is  inversely  proportional  to  turbulence  strength,  becoming  even  smaller  in 
strong  turbulence,  and  the  outer  scale  is  directly  proportional  to  turbulence  strength, 
becoming  larger  in  strong  turbulence. 

It  is  the  random  formation  and  dissipation  of  these  eddies  that  is  responsible 
for  the  stochastic  nature  of  the  atmospheric  IOR.  While  the  mean  value  of  the  IOR 
of  air  no  is  approximately  1,  it  varies  slightly  with  temperature,  pressure,  and  humid¬ 
ity.  These  conditions  are  slightly  different  for  each  eddy.  By  suppressing  the  time 
variations,  the  index  of  refraction  can  be  expressed  as 

n(R)  =  no  +  ni(R)  (2.1) 

where  ni(R)  is  the  deviation  of  the  IOR  from  the  mean  at  some  location  R.  Setting 
no  =  1,  and  writing  n(R)  in  terms  of  temperature  and  pressure  dependence  yields 

n(R)  =  1  +  77.6  x  10~6(1  + 7.52  x  10~3A“2):%^^  (2.2) 

T(R) 

«  1  +  79  x  10-6P-1R)  (2.3) 

T(R)  v  ’ 

where  A  is  the  optical  wavelength  in  microns,  Patm  is  the  pressure  in  millibars  and  T 
is  the  temperature  in  Kelvin.  The  approximation  in  the  second  line  assumes  visible 
wavelengths. 

Statistically,  the  atmospheric  IOR  deviates  from  no  as  a  Gaussian  (or  normal) 
distribution  with  a  slowly  varying  mean.  Though  not  strictly  stationary,  the  atmo¬ 
sphere  does  possess  stationary  increments  which  make  it  well  suited  to  description  by 
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structure  function  Hn(R),  discussed  in  more  detail  later.  Given  that  the  ensemble 
average  of  n(R),  denoted  (n(R)),  is  such  that  (n(R))  =  m(R),  the  spatial  auto¬ 
covariance  function  can  be  expressed  as 

Bn( R)  =  <[n(Ri)  -  m(R1)][n*(R2)  -  m*( R2)]).  (2.4) 

Given  the  further  assumption  of  statistical  homogeneity,  and  substituting  eqn.  (2.1) 
while  noting  that  (ni(R))  =  0,  the  covariance  expression  simplifies  to 

Bn{  R)  =  (ni(Ri)ni(Ri  +  R))  (2.5) 

where  R  =  |R2  —  Ri|.  [1,  Ch.2] 

2.2  Correlated  Random  Variables 

Random  variable  are  correlated  if  they  are  non-orthogonal,  that  is,  if  E[XY]  ^  0 
for  any  random  variables  X  and  Y .  The  covariance  [5,  sec.  4.7]  of  any  two  random 
variables  X  and  Y  is 


COV (X,  Y)  =  E[(X  ~  E[X])(Y  -  E[Y])]  (2.6) 

=  E [XY  -  XE [Y]  -  YE [X]  +  E [X] E [Y]  (2.7) 

=  E[XY]  -  2 E[X]E[Y]  +  E[X]E[Y]  (2.8) 

=  E[XY]  -E[X]E[Y].  (2.9) 


However,  if  X  and  Y  are  independent  (and  thus  uncorrelated),  then  E[XY]  = 
E[X\E[Y]  and  so  COV (X,  Y)  =  0. 

Random  processes,  also  known  as  a  stochastic  processes  have  a  detailed  defi¬ 
nition,  comprehensively  covered  in  [5,  Ch.  6].  They  are  closely  related  to  random 
variables  and  can  be  described  using  the  same  language  of  statistics. 
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The  statistics  of  a  random  process  can  be  functions  of  time,  location,  or  other 
non-random  variables.  Let  t  be  an  independent,  non-random  variable  (possibly  a 
vector),  and  X(t)  be  a  random  process. 

The  mean  of  X(t),  denoted  by  nix,  is 

/OO 

xfX(t){x)dx  (2.10) 

■OO 

where  fx(t)(x )  is  the  probability  distribution  function  (pdf)  of  X(t). 

The  autocorrelation  of  X(t),  denoted  R\,  is 

/OO  /‘OO 

/  xyfX(tl),x{t2)(x,y)dxdy  (2.11) 

oo  J  —  OO 

where  fx[ti),x{t2){x,y)  is  the  second  order  pdf  of  X(t). 

The  autocovariance  of  of  X(t),  denoted  Bx,  can  be  expressed  in  terms  of  the 
mean,  as 

Bx(ti,t2)  =  E[{X(ti)  —  mx(ti)}{X(t2)  —  ‘mx{t2)}]  (2-12) 

and  also  in  terms  of  the  autocorrelation,  as 

Bxituh)  —  Rxitiih)  ~  mx(ti)'mx(t2)-  (2-13) 

Note  that  Bx(ti,t2)  =  Rx(t i,t2)  for  zero  mean  processes. 

For  a  stationary  random  process,  or  one  that  can  be  described  in  stationary 
independent  increments,  the  autocorrelation  can  be  expressed  as  Rx{t)  in  terms  of 
the  offset  r  =  t2  —  t\.  The  Power  Spectral  Density  (PSD)  of  a  random  process  X(t) 
is  defined  as  the  Fourier  Transform  [4,5]  of  the  autocorrelation  Rx(t).  The  PSD, 
denoted  $>x,  then  becomes 

/OO 

Rx(r)e-^Tdr  (2.14) 

■OO 
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Through  the  use  of  Parseval’s  Theorem,  Pave  the  average  power  of  X(t)  can  be 
expressed  as 


Rave  =  E[X2(t )] 

(2.15) 

=  Rx(  0) 

(2.16) 

=  r*xm 

j —00 

(2.17) 

2.3  Structure  Function 

The  structure  function  is  the  expected  variance  (or,  equivalently,  mean  square 
difference)  of  the  random  process  with  respect  to  the  separation  distance.  Mathemat¬ 
ically,  the  structure  function  is  a  difference  function 

-Dn(R-i)  R2)  =  Ti[(n(Ri)  —  n(Ri  +  R))2]  (2-18) 

where  R  =  R2  —  Ri-  In  a  statistically  homogeneous  medium,  the  structure  function 
is  related  to  the  covariance  function  as  follows: 

Ai(R)  =  2[Bn(0)  -  Bn( R)]  (2.19) 

The  additional  assumption  of  isotropy  of  the  propagation  medium  renders  the  statis¬ 
tics  rotationally  invariant  (translational  invariance  is  given  by  homogeneity)  and  al¬ 
lows  the  structure  function  to  be  expressed  purely  as  a  function  of  the  magnitude  of 
the  offset, 

Dn{R)  =  2[Bn{0)  -  Bn(R)}  (2.20) 

where  R  =  |R|. 

There  are  advantages  inherent  in  the  structure  function  that  make  it  preferable 
to  the  covariance  for  this  research.  First,  due  to  the  spectrum  discussed  later,  the 
covariance  has  a  singularity  at  R  =  0  which  the  structure  function  avoids.  Second, 
the  structure  function  is  a  mean  difference  metric,  so  Dn( 0)  will  always  be  zero,  and 


Dn(R  7^  0)  will  be  directly  related  to  the  similarity,  or  correlation,  of  the  values  of 
n(Ri)  to  those  of  n(Ri  +  R),  making  it  an  intuitive  error  metric.  Third,  the  structure 
function  exists  for  non-st  at  ionary  random  process  that  have  stationary  increments, 
whereas  the  covariance  function  is  not  well  defined  for  such. 

The  atmospheric  IOR  structure  function  using  the  Kolmogorov  power  spectrum 
(see  section  2.4)  is  given  in  terms  of  C2  as 


Dn(R)  = 


(2.21) 


C2R2/3,  10CR^:L0 

0%l-4/3R2,  R  <  l0 

where,  as  above,  Z0  is  the  inner  scale  of  the  inertial  subrange.  From  observation, 
Dn  should  increase  exponentially  with  offset  distance,  Dn  oc  i?2/3,  within  the  inertial 
subrange. 


2.4  Power  Spectra 

The  power  spectrum  of  a  process  can  give  enormous  insight  into  the  nature  of 
the  process,  as  it  provides  information  about  the  frequency  content.  This  has  a  direct 
application  in  evaluating  phase  screen  efficacy  because  different  techniques  are  suited 
to  different  ends  of  the  spectrum.  For  instance,  Zernike  polynomials  efficiently  de¬ 
scribe  the  low  frequency  portion  of  the  spectrum  but  can  become  unwieldy  when  the 
upper  frequencies  need  to  be  accurately  modelled;  and  it  is  just  the  opposite  for  the 
FFT  technique.  Because  optical  systems  must  perform  across  the  entire  spectrum, 
phase  screens  must  be  designed  to  include  the  entire  spectrum.  Therefore,  it  is  con¬ 
venient  to  express  the  descriptive  statistics  of  both  phase  screens  and  the  atmosphere 
which  they  simulate  in  terms  of  the  power  spectrum  density  (PSD)  <hn. 

First,  the  power  spectral  density  is  itself  a  Fourier  transform  of  the  autocorrela¬ 
tion  (or  autocovariance  if  the  process  is  zero  mean)  function  of  a  random  process  [1,5] 

*n(K)  =  J  J  f  Bn(R)exp(-jK  •  R )d3R  (2.22) 
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where  k  =  |K|  is  the  scalar  wave  number.  Assuming  that  the  autocorrelation  is 
homogeneous  and  isotropic  the  above  equation  becomes 


(2.23) 


Using  the  inverse  Fourier  Transform  to  express  the  autocovariance  in  terms  of  the 
power  spectrum  yields 


(2.24) 


Using  the  relationship  between  the  structure  function  and  the  autocovariance,  Dn  can 
also  be  expressed  in  terms  of  the  PSD 


(2.25) 


There  are  several  power  spectrum  models,  all  of  which  have  been  derived  em¬ 
pirically  from  observations  of  the  atmosphere  and/or  unit  analysis  of  relevant  quan¬ 
tities  [1]. 

2-4-1  Kolmogorov.  Kolmogorov’s  PSD  is  the  simplest,  but  only  accurate 
within  a  bandwidth  proportional  to  the  inertial  subrange. 


<hn(K)  =  0.033C^k  u/3,  1  /Lq  «k«  1/Z0 


(2.26) 


2-4-2  Tatarski.  Tatarski’s  model  extended  the  accuracy  of  the  model  to 
high  frequencies,  but  is  still  inaccurate  for  k  less  than  1/To- 


<hn(ft)  =  0.033 C2k  u/3exp(-K2 / n2m),  «  >  1/T0 


(2.27) 


Note  that  =  5.92 /l0. 
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10° 


Figure  2.1:  Atmospheric  power  spectra  models  based  upon 

empirical  observations.  Parameters:  /0  =  0.01m,  L0  =  10m 


The  Von  Karmen  PSD  model  is  of  use  across  the  entire 


2-4-3  Von  Karmen. 
spectrum 

*„(«)=  0.033C?K-11/3pt^^,  0  <  k  <  oc 
and  Km  =  5.92/Zq  and  k0  =  1/L0  or  k0  =  2ir/L0. 


(2.28) 


Modified  Atmospheric.  The  modified  Atmospheric  spectrum  is  the 
most  complicated  of  the  PSD  models  discussed  herein.  It  largely  follows  the  Von 
Karmen  PSD,  but  allows  for  the  ‘bump’  observed  around  1/Iq,  as  follows 

$„(«)  =  0.033C^[1  +  1.802(ac/acz)  -  0.254(/c//ct)7/6]  ,  0  <  k  <  oc  (2.29) 

where  /y  =  3.3/Zq  and  kq  =  1/To  or  k0  =  2tt/L0. 
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2.5  Atmospheric  Turbulence  Descriptors 

Optical  atmospheric  models  are  the  subject  of  ongoing  research.  Most  models 
include  the  concept  of  the  atmospheric  index  of  refraction  structure  constant 
denoting  the  intensity  of  atmospheric  turbulence.  The  integrated  C%  profile  is  defined 
as 

cl  =  J  Cl(z)dz  (2.30) 

where  2  is  the  optical  path.  can  change  drastically  with  time  of  day,  location  and 
weather  conditions.  However,  for  relatively  clear  conditions,  an  reasonable  C\  value 
is  on  the  order  of  1CF16. 

For  a  plane  wave,  still  using  the  Kolmogorov  power  spectrum,  the  phase  struc¬ 
ture  function  D^R)  is  related  to  the  integrated  as 

D^R)  =  2.91  k2R5/3C*  (2.31) 

where  k  is  the  wave  number.  The  atmospheric  coherence  radius  po  is  defined  as 
the  point  when  the  D^R)  =  2.  This  corresponds  to  the  separation  at  which  the 
autocorrelation  is  1/e.  The  Fried  radius  r0  is  related  to  po  as 


t  0  —  2.1po  (2.32) 

The  Fried  radius  is  used  together  with  the  diameter  of  a  system  optic  to  provide  a 
single  widely  used  system  descriptor,  D/tq. 

In  terms  of  C/,  ro  can  be  derived  by  noting  the  relationship  between  the  phase 
and  atmospheric  structure  functions,  in  a  homogeneous,  isotropic  medium,  is 

D^R)  =  2.91  k2RDn(R)  (2.33) 
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Figure  2.2:  Theoretical  phase  structure  functions  for  various 
D/r0  values. 


Using  the  equation  above  along  with  equation  2.21,  r0  of  a  plane  wave  can  be  expressed 
as 


r0 


0.185 


47T2 

k2C2 


3/5 


(2.34) 


This  in  turn  allows  the  D $  to  be  expressed  as 


D^R)  =  6.88(i?/r0)5/3 


(2.35) 


where  the  proportionality  D $  to  i?5/3  is  particularly  evident. 

A  full  derivation  of  the  phase  and  atmospheric  IOR  structure  functions  can  be 
found  in  Roggemann  and  Welsh’s  Imaging  Through  Turbulence  [9,  sec.  3. 3, 3. 4], 


2. 6  Types  of  Phase  Screens 

Phase  screens  are  used  in  wave  optics  models  to  simulate  the  correlated  ran¬ 
dom  phase  change  that  occurs  during  propagation  of  light  through  turbulent  media, 
such  as  the  atmosphere.  There  are  three  methods  commonly  used  to  generate  phase 
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Table  2.1:  Zernike  Polynomials. 


n 

m 

i 

Polynomial 

Name 

0 

0 

1 

1 

Piston  (Ignored  herein) 

1 

1 

2 

2  r  cos  9 

Tip 

1 

1 

3 

2  r  sin  9 

Tilt 

2 

0 

4 

3.464r2  -  1.732 

Defocus 

2 

2 

5 

2.449r2  sin  29 

Astigmatism  1 

2 

2 

6 

2.449r2  cos  29 

Astigmatism  2 

3 

1 

7 

(8.485r3  -  5.657r)  sin  9 

Coma  1 

3 

1 

8 

(8.485r3  -  5.657r)  cos  9 

Coma  2 

3 

3 

9 

2.828 r3  sin  3# 

3 

3 

10 

2.828r3  cos  3  9 

4 

0 

11 

3.416r4  -  13.416r2  +  2.236 

Spherical  Abberation 

screens,  Zernike  polynomials,  Fast  Fourier  Transforms  (FFT),  and  Fourier  Series  (FS). 
There  is  also  a  modification  of  the  FFT  technique  called  Sub-Harmonic  Frequency 
Expansion  (SHFE)  which  is  worth  mentioning,  but  is  not  covered  in  depth  due  to  its 
similarity  to  both  the  FFT  and  the  FS  methods.  Each  of  these  phase  screen  genera¬ 
tion  methods  relies  upon  random  magnitudes  applied  to  sets  of  basis  functions.  For 
the  Zernike  method,  the  functions  are  a  set  of  indexed  polynomials,  while  the  FFT 
and  FS  methods  rely  upon  sinusoids. 


2.1  Zernike  Polynomials 

Zernike  phase  screens  are  generated  as  a  weighted  set  of  polynomials  [9].  Each 
Zernike  polynomial  is  continuous  and  describes  a  unique  pattern  of  abberation  as 
may  be  found  in  a  lens  of  unit  diameter.  The  set  Zernike  polynomials  form  a  normal, 
orthogonal  basis  set  thus,  for  any  Zernike  polynomials  and  Zj 

27r  f1  [  0  if  i  7^  j 

/  Zi(r,6)Zj(r,6)drd6  =  <  (2.36) 

Jo  I  1  if  i  —  j 
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The  equations  in  table  2.1  are  derived  from  the  more  general  equations 


teven(Xi 

=  y/2(n  +  l)R™(r)  cos(m6)),  m  ^  0 

(2.37) 

i^oddix 5 

=  a/2 (n  +  1  )R™(r)  sin (m9),m  ^  0 

(2.38) 

Zi{r,e) 

=  R°n,m  =  0 

(2.39) 

where  mn  and  n  are  the  indices  of  azimuthal  order  and  radial  order,  respectively,  and 
r  <  1.  If  a  screen  of  other  than  unit  radius  (r  >  1)  is  required,  the  following  transform 
is  used:  For  a  screen  of  desired  radius  R ,  with  absolute  position  as  r,  and  normalized 
position  p  =  r/R,  the  phase  can  be  expressed  as 


OO 


4>(Rp, 

=  ^2aiZi(p,d) 

2—1 

(2.40) 

00",/) 

OO 

=  ^  ciiZi  ,  o'j 

2—1 

(2.41) 

Note  that  table  2.1  also  uses  the  index  i,  which  is  used  in  applying  the  Zernike 
coefficient  Oj.  The  subscript  i  is  unique  to  an  individual  polynomial  Zt  (as  n  and  m 
are  not)  and  allows  one  to  specify  the  power  of  each  separate  abberation  within  the 
system  once  the  polynomials  have  been  generated. 

The  Zernike  coefficients  a*  relating  to  to  Zernike  polynomials  Z%  can  be  recovered 
from  a  given  phase  screen  due  to  their  orthonormality,  using  the  following  equation 

=  fW(r)0(r)Zj(r)dr 

f  W (r) | Zi (r) \2dr  1  '  ’ 

where  W (r)  is  a  windowing  function  that  defines  the  extent  of  the  Zernike  polynomials, 
( f>(r)  is  the  zero  mean  input  phase  map  and 


/ 


W(r)dr  = 


1 
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Figure  2.3:  Zernike  Polynomials  2  through  10.  Note  that  the 
Zenike  basis  set  assumes  a  unit  circle  and  must  be  scaled  if  that 
is  not  the  case. 

The  assumption  that  <fi(r )  is  zero  mean  implies  that  piston  error  (Zi)  is  ignored, 
which  is  practical  as  it  corresponds  to  a  non-measurable  constant  time  delay  from  the 
source  to  the  detector.  If  the  screen  is  not  zero  mean,  it  must  be  normalized  before 
equation  2.42  is  applied.  cq  can  be  recovered  by  noting  that  it  is  equal  to  the  mean 
of  the  phase  before  normalization. 

The  total  aberrant  power  (i.e.  the  average  mean  square  wavefront  error)  of  a 
Zernike  phase  screen  can  be  expressed  by  the  sum  of  its  squared  Zernike  coefficients, 
excluding  a±  for  the  reasons  given  above. 


=  J  W{r)(j)2{r)dr  (2.43) 

OO 

=  I2-44) 

2 
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Table  2.2:  Zernike  Covariance  Ma¬ 

trix  (Empty  squares  represent  zeros). 


2 

3 

4 

5 

6 

7 

8 

9 

10 

2 

0.448 

-0.0141 

3 

0.448 

-0.0141 

4 

0.0232 

5 

0.0232 

6 

0.0232 

7 

-0.0141 

0.00618 

8 

-0.0141 

0.00618 

9 

0.00618 

10 

0.00618 

Table  2.2  shows  the  Zernike  covariance  matrix  for  a  system  D/tq  =  1,  denoted 
A ij,  using  the  Komogorov  spectrum.  The  covariance  matrix  for  any  given  D/ro  value 

is 


Bz.Zj  =  A  ^(D/ro)5'3  (2.45) 

Scaling  the  covariance  matrix  and  applying  it  to  the  set  of  polynomials  are  key 
steps  in  construction  a  Zernike  phase  screen,  as  is  further  discussed  in  chapter  3. 

The  Zernike  method  builds  phase  screens  from  low  frequency  to  high,  so  the  first 
few  polynomial  describe  the  abberations  of  greatest  power.  As  noted  by  Roggeman 
and  Welsh  [9]  nearly  80%  of  the  aberrative  power  of  an  optical  system  is  contained 
within  the  tip  and  tilt  ( and  Z3)  error  alone.  This  is  advantageous  for  ‘ball-park’ 
simulations,  when  an  adequate  phase  screens  may  be  constructed  with  only  a  few 
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polynomials.  The  disadvantage  arises  when  attempting  more  precise  high-frequency 
models,  as  they  can  require  a  large  number  of  polynomials  to  be  generated. 

2.8  Fourier  Phase  Screens 

The  discrete  Fourier  transform  is  the  representation  of  a  signal  as  the  sum  of  its 
frequency  content  [4,6,8].  The  result  is  a  complex  valued  function  of  the  frequency 
variable  k  =  (kx,ky): 


F(kx,  ky ) 


TV— 1  TV-1 

j2J2f^y^mxkx+vky) 

£=0  y= 0 

(2.46) 

F r{]^,x'i  ky)  F  j  F i(kx,  ky ) 

(2.47) 

\F(kx,ky)\ejZF^kv) 

(2.48) 

Both  the  FFT  and  the  FS  methods  use  this  transform,  as  does  SHFE.  However, 
the  difference  between  the  methods  is  readily  discussed  using  the  inverse  transform. 
The  two-dimensional  inverse  transform  represents  a  signal  f(x,y )  as  the  sum  of  its 
frequency  content  F(kx,  ky) 

f{x,  y)  =  Y,Y, F (**’  ky)e~^xk°+^  (2.49) 

kx  f^y 

The  FFT  method  uses  frequencies  kx  and  ky  at  regular  intervals  on  a  cartesian  grid, 
such  that  dkx  =  dky  for  all  kx  and  ky  (see  figure  2.4).  The  FS  method  presented 
here  uses  frequencies  at  logarithmic  intervals  on  a  polar  grid  (figure  2.5).  The  SHFE 
method  uses  the  same  base  frequency  grid  as  the  FFT,  but  iteratively  expands  the 
lower  frequencies  (figure  2.6)  to  build  more  low  frequency  content. 

2. 8. 1  Fast  Fourier  Transform.  By  evenly  spacing  the  component  frequencies 
on  a  cartesian  grid  and  using  code  that  has  been  extensively  optimized,  the  FFT 
method  requires  only  N  log  N  operations,  reducing  the  number  of  operations  required 


18 


0 


Figure  2.4: 


Figure  2.5: 


o 


The  base  frequency  grid  of  an  FFT  phase  screen. 


x  io4 


The  base  frequency  grid  of  an  FS  phase  screen. 
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Figure  2.6:  The  zero  frequency  pixel  of  an  FFT  phase  screen 
after  three  iterations  of  Sub-Harmonic  Frequency  Expansion. 

by  an  order  of  N/logN.  The  atmospheric  phase  is  modelled  via  a  power  spectrum 
$(k)  used  as  a  filter.  The  filter  is  applied  to  Gaussian  random  variables  to  create  the 
weighting  function  for  each  frequency  component,  and  the  Fourier  transform  is  taken. 

Sfft( f)  =  F-'iGWy/Qtf)}  (2.50) 

where  G(k)  is  a  matrix  of  unit  variance  Gaussian  random  variables  and  is  the 
inverse  Fourier  transform  operator.  If  G(k)  is  complex  Gaussian, 

Sfft(t)  =  SR(r)  +  jSj(r) 

where  Sr  and  Sj  are  independent.  This  yields  two  distinct  phase  screens  0i(r)  = 
SR{f)  and  (j)2{f)  = 

There  are  two  negative  characteristics  associated  with  the  FFT  frequency  grid. 
First,  all  frequencies  lower  than  the  sampling  frequency  fs  are  lumped  into  the  DC 
component  (see  figure  2.4.  Due  to  the  fact  that  atmospheric  turbulence  power  is 
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Figure  2.7:  The  nine  lowest  frequencies  of  an  FFT  phase 

screen.  The  diagonal  frequency  points  are  farther  from  zero 
than  the  on- axis  points  by  a  factor  of  y/2. 

highest  for  low  frequencies,  this  may  lower  fidelity.  Second,  the  frequency  steps  are 
larger  on  the  diagonal  than  on  either  of  the  axes  by  a  factor  of  y/2  (see  figures  2.7 
and  2.8). 

2.8.2  Fourier  Series.  The  FS  phase  screen  generated  on  a  logarithmically 
scaled  polar  grid,  rather  than  an  equidistant  rectangular  grid.  The  logarithmic  scaling 
allows  for  a  higher  density  of  frequencies  near  zero  while  not  requiring  the  same  density 
at  the  higher  frequencies  (figure  2.9).  Figure  2.10  is  the  magnified  center  of  figure  2.5, 
showing  the  low  frequency  components.  No  that  the  low  frequency  grid  is  divided  into 
multiple  low  frequencies.  The  polar  distributions  ensures  that  the  frequency  step  will 
be  the  same  regardless  of  direction;  eliminating  the  y/2  difference  seen  on  the  FFT 
diagonals. 

2.8.3  Sub-Harmonic  Frequency  Expansion.  In  order  to  avoid  generating 
excess  high  frequency  data  while  retaining  the  speed  of  the  FFT,  the  Sub-Harmonic 
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Figure  2.8:  Example  of  the  nine  lowest  frequencies  within  an 
FFT  phase  screen,  corresponding  to  figure  2.7.  The  effect  of  the 
v/2  factor  in  the  frequency  step  of  the  diagonal  elements  can  be 
seen  in  the  closer  spacing  of  the  lines  in  (a),(c),(g)  and  (f). 


Figure  2.9:  Frequency  domain  locations  of  the  sinusoids  con¬ 
tained  within  the  FS  phase  screen  on  log-log  and  semi-log  axes. 
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Figure  2.10:  Frequency  domain  locations  of  the  sinusoids  con¬ 
tained  within  the  FS  phase  screen. 

Frequency  Expansion  (SHFE)  can  be  used  to  expand  only  the  low  frequency  terms  of 
an  FFT  phase  screen.  This  is  done  by  dividing  the  zero  frequency  pixel  of  the  FFT 
screen  into  increasingly  finer  grids  (2.11).  Each  new  frequency  is  used  to  generate  a 
two-dimensional  sinusoid  which  is  then  weighted  with  a  random  added  to  the  original 
FFT  phase  screen.  Due  to  the  decrease  in  frequency,  the  phase  screen  will  be  to  small 
to  contain  an  entire  sinusoidal  period,  causing  the  screen  to  gain  a  non-zero  mean. 
Therefore,  the  phase  screen  may  need  to  be  re-normalized  after  SHFE. 

The  SHFE  method  combines  the  optimized  code  of  the  FFT  method  with  some 
of  the  individual  frequency  specihcation  of  the  FS  method.  It  is  more  cumbersome 
than  the  pure  FFT  method  and  still  retains  some  of  the  flaws,  such  as  the  a/2  difference 
between  axial  and  diagonal  frequency  steps  (figure  2.12).  The  SHFE  frequency  spread 
is  not  as  efficient  as  that  of  the  FS  method,  even  though  it  does  increase  frequency 
density  near  zero.  The  SHFE  method  can  be  viewed  as  a  hybrid  combination  of  FFT 
and  FS  methods.  Therefore,  though  worth  mentioning  for  completeness,  SHFE  will 
not  be  evaluated  individually  herein. 
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Figure  2.11:  The  low  pixels  of  an  FFT  frequency  grid  after 
one  iteration  of  SHFE.  The  expanded  frequencies  are  denoted 
with  an  asterisk. 


Figure  2.12:  Phase  effects  examples  corresponding  to  the  ex¬ 
panded  frequencies  from  figure  2.11.  Note  that  they  are  lower 
frequency  than  those  of  2.8  but  still  contain  the  y/2  difference 
on  the  diagonals. 
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III.  Methods  of  Comparison 


Three  types  of  phase  screens  are  generated  -  Zernike,  FFT,  and  FS  -  in  various 
sizes.  Each  of  which  must  then  evaluated  for  time  efficiency  and  correspondence 
to  theory.  This  chapter  outlines  the  evaluation  methods  used  herein.  Section  3.1 
covers  the  methods  used  to  generate  all  three,  in  the  above  order.  Section  3.2  covers 
what  needs  to  be  considered  when  deciding  the  size  of  phase  screen  to  generate. 
Section  3.3  discusses  how  the  generation  time  of  each  method  is  analyzed.  Section  3.4 
contains  information  about  the  frequency  content  of  the  Zernike  polynomials.  Section 
3.5  is  a  brief  summary  of  the  use  of  the  structure  function. 

3.1  Generating  Phase  Screens 

3.1.1  Zernike.  A  Zernike  phase  screen  is  generated  by  applying  weighted 
coefficients  to  the  set  of  Zernike  polynomials.  The  coefficient  values  are  determined  by 
the  Zernike  covariance  matrix  for  the  given  D/r o.  The  random  samples  are  calculated 
by  using  the  cholescky  factorization  via  the  MATLAB®  ‘choh  command  [6] 

A  =  CHOL(Kay  *  G  (3.1) 

where  Ka  is  the  Zernike  covariance  matrix,  G  is  a  vector  of  unit  variance,  zero  mean 
Gaussian  random  variables  and  A  is  the  resulting  vector  of  random,  correlated  Zernike 
coefficients.  Note  that  the  Zernike  covariance  matrix  for  any  D/tq  can  be  generated 
by  multiplying  the  covariance  matrix  for  D/r0  =  1  by  (D /r0)^5^3\  This  allows  one 
Zernike  covariance  matrix  to  be  loaded  from  memory  and  used  with  only  a  multi¬ 
plicative  change. 

The  MATLAB®  routine  MakeZernikePhzScrn.m  generates  an  N  x  N  pixel,  zero 
mean  phase  screen  with  a  covariance  corresponding  to  a  given  D/tq.  ft  takes  as  input 
four  parameters:  the  system  D/tq  value,  the  size  (in  pixels)  of  the  screen,  the  square 
Zernike  covariance  matrix  for  D/ro  =  1  (Kao),  and  the  number  of  Zernike  polynomials 
to  be  used.  The  process  of  generating  a  Zernike  phase  screen  can  be  summarized  as 
follows: 
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Figure  3.1:  A  Zernike  Phase  Screen  constructed  using  Zernike 
polynomials  2  through  81.  N  —  100  pixels  per  side. 

•  Calculate  given  Zernike  polynomials  Z*  on  a  unit  circle. 

•  Load  Zernike  covariance  matrix  Ka  and  normalize  it  to  D/tq 

•  Calculate  random  weights  a*  and  apply  them  to  Zr. 

•  Sum  Zi  s,  after  any  necessary  scaling. 

There  are  three  items  of  particular  note  when  generating  a  Zernike  phase  screen 
using  MakeZernikePhzScrn.m.  First,  the  dimensions  of  Ka0  must  be  at  least  as  large 
as  the  number  of  Zernike  polynomials  to  be  used,  as  it  does  in  fact  describe  their 
covariance.  Second,  the  first  degree  Zernike  polynomial  (describing  time  delay,  or 
‘piston’  error)  is  not  present  in  the  covariance  matrix  and  not  used  in  constructing 
the  phase  screen.  The  user  can  compensate  for  this  by  adding  a  constant  to  the  screen, 
making  it  a  non-zero  mean  matrix.  Third,  as  can  be  seen  in  figure  3.1,  Zernike  phase 
screens  are  circular.  Therefore  the  usable  screens  size  is  only  7rfV2 /4  rather  than  N2 
which  is  the  size  of  the  matrix  generated. 
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3.1.2  FFT.  The  MATLAB®  code  Ma.ke_FFTPhzScrn.m  generates  a  phase 
screen  corresponding  to  a  particular  atmospheric  spectral  model  $  at  a  given  C2. 
The  input  spectra  $  must  be  the  same  dimensions  as  the  desired  phase  screen.  A 
peculiarity  of  the  FFT  method  is  the  continuity  across  edges.  This  leads  to  higher  than 
theory  correlation  of  values  on  opposite  sides  of  the  screen,  causing  non-physical  edge 
effects  (Figure  3.2).  This  is  an  inherent  result  of  the  Fourier  transform  method  and 
can  be  compensated  for  by  generating  a  larger  than  necessary  screen  and  discarding 
the  edges.  In  summary,  to  build  an  FFT  phase  screen 


•  Build  a  two-dimensional  cartesian  frequency  grid  (shown  here  using  Von  Karmen 
PSD) 


K  = 


-11/12 

exp  [-0.56 (l0/D)2(kl  +  k2y)\ 


•  Generate  random  weights  G  as  a  set  of  zero  mean,  unit  variance,  complex 
Gaussian  random  numbers. 


Take  the  Fourier  transform  of  the  weighted  frequency  grid  kG  and  multiply  by 


0.1517  (  — 

r  0 


(5/6) 


Take  the  real  and  imaginary  parts  of  the  result  to  form  two  distinct  phase 


screens. 


3.1.3  FS.  An  FS  phase  screen  is  generated  by  first  constructing  a  kernel  of 
two-dimensional  sinusoids  with  frequencies  relating  to  logarithmically  spaced  points 
on  a  polar  grid.  Using  the  multiplicative  spacing  factor  Q,  the  number  of  frequencies 
Nk  is  inversely  proportional  to  the  value  of  Q,  such  that 


NK  =  Ne 


log 


14.2Uq 

trl0 


log  Q 


(3.2) 
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Figure  3.2:  An  FFT  Phase  Screen  constructed  using  Von  Kar¬ 
men  frequency  spectrum. 

where  Ng  is  the  number  of  equal  angles  in  each  polar  grid.  Using  a  triply  nested 
for  loop,  and  the  mean  PSD,  the  value  of  each  sinusoid  at  each  pixel  is  calculated, 
yielding  a  three-dimensional  kernel  size  of  Sk  =  Nx  x  Ny  x  NK.  Generating  square 
phase  screens,  as  is  done  here,  N  =  Nx  =  Ny,  so  the  kernel  size  is 

Sk  =  N2x  Nk  (3.3) 

Randomness  is  applied  to  the  kernel  by  multiplying  by  a  zero  mean,  unit  vari¬ 
ance,  complex  Gaussian  random  number.  The  real  and  imaginary  parts  are  taken  to 
form  two  sets  of  data,  then  each  set  is  summed  and  the  mean  is  subtracted,  forming 
two  distinct  random  phase  screens.  In  summary,  to  form  an  FS  phase  screen 

•  Generate  a  kernel  sinusoids  with  frequencies  relating  to  a  logarithmic  polar  grid, 
using  the  mean  PSD. 

•  Weight  each  sinusoid  with  a  zero  mean,  unit  variance,  complex  Gaussian  random 
number. 
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Figure  3.3:  An  FS  Phase  Screen  constructed  using  a  logarith¬ 
mic  kernel. 

•  Use  the  real  and  imaginary  parts  of  the  result  to  for  two  distinct  phase  screens. 

Note  that  the  itemized  steps  for  FS  greatly  resemble  those  of  the  FFT.  The  main 
difference  is  that,  due  to  the  frequency  grid  used,  the  FS  routine  cannot  make  use  of 
the  optimized  code  in  MATLAB®  and  must  instead  construct  each  sinusoid  a 

pixel  at  a  time.  The  result  (figure  3.3)  contains  a  greater  spectrum  representation  of 
low  frequency  content  than  in  typical  FFT  techniques.  It  also  avoids  the  continuity 
across  edges  that  causes  some  non-physical  edge  effects. 

3. 2  Sizing  Considerations 

When  generating  phase  screens  as  part  of  a  model  of  physical  phenomena  the 
physical  area  represented  is  important.  The  physical  size  is  specified  by  the  sample 
period  Ax  and  the  number  of  pixels  across  N,  yielding  a  total  area  A  of 

A  =  NAx  x  NAx  =  N2Ax2  (3.4) 
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As  mentioned  above,  both  Zernike  and  FFT  phase  screens  have  a  useable  size 
somewhat  less  than  the  N  x  N  pixel  grid  generated.  For  FFT  phase  screens,  this  is 
because  the  method  itself  causes  the  opposite  edges  of  the  generated  phase  screen  to 
be  highly  correlated,  which  is  contrary  to  observation  of  the  atmosphere.  This  is  a 
problem  that  can  be  overcome  by  generating  a  phase  screen  large  than  what  is  needed 
and  using  only  the  center  portion  for  simulation. 

In  the  case  of  Zernike  phases  screens,  the  size  discrepancy  is  the  result  of  the 
definition  of  Zernike  polynomials.  Each  Zernike  polynomial  is  defined  only  on  the  unit 
circle,  and  must  be  scaled  to  represent  a  larger  area  using  the  following  transform 

r  =  Rap  (3.5) 

where  r  is  the  absolute  position,  Ra  is  the  constant  radius  of  aperture  and  p  is  the 
normalized  position.  However,  even  once  scaled,  the  corners  of  the  square  screen  will 
still  be  zero,  leaving  a  represented  area  of 

Az  =  j N2Ax 2  (3.6) 

FS  screens  avoid  both  limitations,  being  defined  as  sinusoids  upon  an  infinite 
grid  and  having  a  distribution  of  frequency  components  that  does  not  artificially 
correlate  the  outer  edges  of  the  finite  screen. 

3.3  Generation  Time  Analysis 

Each  phase  screen  has  some  parameter  (s)  that  may  be  generated  once  at  the 
beginning  of  given  scenario  and  used  to  make  any  number  of  screens.  This  initial 
generation  time  T*  can  be  separated  from  the  individual  screen  generation  time  Ts, 
allowing  the  total  time  used  to  generate  Ns  screens  (Tt)  to  be  expressed  as 

Tt=Ti  +  NSTS.  (3.7) 
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Optimization  of  Tt  must  consider  the  optimal  number  of  phase  screens  needed 
to  effectively  model  a  given  scenario.  As  the  number  of  screens  needed  goes  up,  the 
effect  of  Tt  becomes  less  important.  Also,  it  may  be  feasible  to  complete  all  initial 
calculations  prior  to  the  implementation  the  simulation,  again  lessening  the  emphasis 
upon  initialization  time.  However,  if  one  must  change  conditions  drastically  during 
the  simulation,  T,;  can  have  a  great  impact  upon  the  efficiency  of  the  model. 

3.3.1  Initialization  Time.  In  MATLAB®,  execution  time  can  be  measured 
using  the  tic  and  toe  commands.  In  the  case  of  the  Zernike  PS,  tic  is  used  before 
loading  the  covariance  matrix  and  Zernike  polynomial  index;  toe  is  used  after  they 
are  loaded.  The  Zernike  PS  initialization  time  T^_z  is  then  set  equal  to  toe. 

Similarly,  tic  and  toe  can  be  used  to  record  the  FFT  PS  initialization  time 
Ti-FFT,  and  the  FS  PS  initialization  time  T^ps-  Ti-fft  encompasses  calculating 
the  analytical  spectra  (<f>).  T^fs  encompasses  defining  a  frequency  grid  spacing  and 
generating  a  polar  kernel  of  sinusoids,  one  for  each  frequency  grid  point. 

3.3.2  Per  Screen  Time.  The  per  screen  generation  time  Ts  is  calculated  by 
using  tic  and  toe  to  record  the  time  it  takes  to  construct  N  phase  screens  after  the 
initial  constructs  have  been  formed  (T)),  then  dividing  by  N: 

Ts  =  C Tt  -  Ti)/N 

However,  in  practice  T%  and  Ts  are  recorded  independently  within  the  code. 

3-4  Zernike  Frequency  Analysis 

Taking  the  FFT’s  of  the  Zernike  equations  (2.37)  both  individually  and  as  a 
screen  made  of  of  equally  weighted  polynomials  allows  one  to  compare  the  frequency 
content  of  the  Zernike  phase  screens  with  that  of  the  FFT  and  FS  phase  screens.  This 
is  done  by  generating  an  set  of  Zernike  polynomials,  Z%  through  27,-,  using  non-random 
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coefficients 


a*  =  dj  for  all  i,j 


The  polynomials  are  then  Fourier  transformed,  and  their  power  spectrum  evaluated, 
both  in  one  dimension  and  two.  For  the  purpose  of  cross-sectional  comparisons,  a 
selection  of  polynomials  with  increasing  radial  order  n  is  used.  Results  are  presented 
in  chapter  4. 

3.5  Structure  Function  Analysis 

Being  a  mean  square  difference  function,  the  structure  function  of  a  phase  screen 
can  be  found  as  follows:  For  any  phase  screen  0,  let  <f>{B)  be  the  original  phase  screen 
0(0)  shifted  in  some  direction  by  R  pixels,  then 


D*(R)  =  S[(0(O)  -  ^(fi))2] 


(3.8) 


Implementing  this  for  a  phase  screen  with  values  distributed  on  a  discrete  carte¬ 
sian  grid  means  that  D^R)  is  affected  by  the  direction  of  the  frequency  shift.  If  the 
shift  is  along  either  of  the  axes,  each  step  A R  will  be  a  single  unit  of  separation. 
However,  if  it  is  off  axis,  such  as  along  the  diagonal,  each  step  will  be  greater  than  a 
single  unit 


A  R  =  x/(Ai  +  Ai/)2 


(3.9) 


and  D^(R)  must  be  scaled  accordingly  when  plotted. 
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IV.  Results 


This  chapter  contains  the  results  from  the  methods  described  in  chapter  3.  Sec¬ 
tion  4.1  covers  the  generation  times  for  all  three  methods,  as  well  as  how  the  the 
Zernike  method  and  FS  method  times  are  affected  by  varying  Nz  and  Q,  respectively. 
Section  4.2  addresses  how  the  frequency  content  of  Zernike  phase  screens  compares  to 
that  of  FFT  and  FS  phase  screens.  Section  4.3  covers  how  the  structure  functions  of 
the  various  screens  relate  to  the  theoretical  atmospheric  phase  structure  function.  Fi¬ 
nally,  section  4.3  discusses  quality /usability  comparisons  between  FFT  and  FS  phase 
screens. 

4-1  Time  Results  and  Analysis 

This  section  is  broken  down  into  four  subsections:  the  first  two,  initialization 
time  and  per  screen  time,  give  the  usability  of  each  method  in  terms  of  how  long 
each  one  takes.  The  second  two,  Zernike  time  explanation  and  FS  time  explanation, 
are  included  because  there  is  much  optimization  that  the  routines  lack.  The  FFT 
code  has  already  been  extensively  optimized  because  it  is  used  so  ubiquitously  within 
signal  processing,  and  so  there  is  no  particular  need  to  focus  upon  it.  The  results 
presented  here  are  based  upon  the  performance  of  a  particular  computer,  so  the  times 
will  change  given  different  hardware.  However,  they  should  be  accurate  as  a  relative 
measure  of  performance. 

4-1.1  Initialization  time.  The  initialization  time  T*  includes  different  ele¬ 
ments  for  each  method.  Initialization  of  an  N  x  N  phase  screen  requires  the  following: 
Using  FFT’s,  one  must  initially  calculate  the  NxN  atmospheric  spectra  using  a  model 
such  as  those  discussed  in  section  2.4.  Using  FS,  one  must  generate  the  polar  kernel, 
choosing  each  included  frequency  on  a  two-dimensional  grid.  Each  frequency  point 
corresponds  to  an  N  x  N  sinusoid  in  the  time  domain.  Using  Zernikes,  one  must  load 
both  a  covariance  matrix  and  a  polynomial  index,  then  each  indexed  polynomial  must 
be  generated  in  a  N  x  N  matrix. 
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Figure  4.1:  Initialization  times  for  all  three  types  of  phase 

screens  with  respect  to  the  number  of  pixels  per  screen  edge  N. 

Figure  4.1  shows  the  initialization  times  for  all  three  methods.  Times  were 
calculated  using  84  Zernike  polynomials  and  a  Q  value  of  2.0.  Recall  that  T^fft, 
Ti_FSi  Ti-z  are  the  initialization  times  for  the  FFT,  FS,  and  Zernike  methods,  re¬ 
spectively,  and  Ts_fft ,  Ts_Fs,  Ts_z  are  their  respective  per-screen  times.  Tt_FFT  is 
by  far  the  smallest  and  is  nearly  linear  in  the  range  given.  This  is  due  to  the  small 
amount  of  preparation  needed  for  the  FFT  method,  most  of  which  makes  use  of  opti¬ 
mized  MATLAB®  vector  operations.  T^z  and  T^ps  are  more  interesting.  Both  are 
quadratically  proportional  to  the  number  of  pixels  per  screen  edge  N,  making  them 
linearly  proportionate  to  the  number  of  pixels  per  screen,  N2. 

Ti-ps  is  the  highest  due  to  the  kernel  size  and  the  nested  for  loop  used  in 
generating  it.  The  FS  kernel  size  Sp  is  proportional  to  N2  and  the  for  loop  execution 
time  is  directly  proportional  to  the  kernel  size.  The  FS  initialization  encompasses 
the  generation  of  416  N  x  N  sinusoids.  By  comparison,  T^_z  is  lower  because  it 
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Figure  4.2:  (a)  The  time  to  generate  each  phase  screen  af¬ 

ter  the  initial  calculations  have  been  made,  with  respect  to  the 
square  root  of  the  number  of  pixels  in  the  screen. 

encompasses  the  calculation  of  only  84  N  x  N  polynomials.  T^fft  is  lowest  of  all 
because  it  encompasses  the  calculation  of  only  4  N  x  N  power  spectra. 

Judged  purely  upon  initialization  time,  the  FFT  method  is  by  far  the  most 
efficient.  The  Zernike  method  comes  second  and  the  FS  method  is  the  least  efficient. 

4-1.2  Per-Screen  Time.  The  per-screen  time  Ts  (figure  4.2)  is  taken  from 
those  calculations  that  cannot  be  made  in  advance  for  a  given  set  of  phase  screens. 
For  each  method,  this  involves  the  generation  of  complex  random  weights.  In  the 
case  of  the  FFT,  it  also  involves  taking  the  two-dimensional  fast  Fourier  transform  of 
the  N  x  N  data  matrix  and  separating  the  result  into  its  real  and  imaginary  parts. 
In  MATLAB®,  this  process  has  been  extensively  optimized  using  vector  operations, 
with  excellent  results. 

As  can  be  seen  in  figure  4.2,  Ts_fft  is  the  lowest  of  the  three  per-screen  gen¬ 
eration  times  shown.  Figure  4.3  shows  this  more  clearly  by  ignoring  the  far  larger 
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Figure  4.3:  (a)  The  time  to  generate  each  Zernike  and  FFT 

phase  screen  after  the  initial  calculations  have  been  made,  with 
respect  to  the  square  root  of  the  number  of  pixels  in  the  screen. 

Ts_ps  and  showing  only  Ts_fft  and  TS_Z-  Though  not  entirely  linear  and  increasing 
monotonically,  Ts_fft  follows  the  expected  upward  trend  proportional  to  the  number 
of  operations  required,  iVlogiV.  It  does  not  rise  above  a  tenth  of  a  second  until  after 
N  =  320.  Figure  4.3  also  shows  that  Ts_z  passes  the  tenth  of  a  second  mark  shortly 
after  N  =  128,  and  figure  4.2  shows  that  Ts_fs  rises  above  0.01  s  for  IV  =  100. 

The  irregular  progression  of  Ts_fs  in  figure  4.2  is  particularly  noteworthy.  Al¬ 
though  it  possesses  a  general  quadratically  increasing  trend,  it  is  not  monotonic.  De¬ 
creases  are  evident  from  N  =  128  to  N  =  192  and  again  from  N  =  256  to  N  =  320. 
These  are  most  likely  a  result  of  the  manner  in  which  MATLAB®  allocates  memory: 
one  segment  at  a  time. 

By  default,  MATLAB®  uses  64-bit  (8-byte)  double  precision  floating  point  num¬ 
bers.  The  memory  requirements  M  for  a  kernel  can  therefore  be  calculated  in  kilobytes 
(KB) 

M  =  8SV/1024  =  81Y2]VK/1024  (4.1) 
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using  the  binary  definition  of  1  KB  =  1024  bytes  and  1  MB  =  1024  KB.  For  N  =  100 
and  NK  =  416  (i.e.,Q  =  2.0),  the  memory  required  is  M  =  32,  500  KB  rs  33  MB.  For 
N  =  384  at  the  same  Q  value,  M  =  479232  KB  =  468  MB. 

Complex  numbers  are  treated  as  two  separate  floating  point  numbers,  the  real 
part  ant  the  imaginary  part.  Thus  one  complex  number  actually  takes  128  bits  (or  16 
bytes)  of  memory.  Therefore,  when  the  random  weights  are  applied  to  the  kernel  to 
generate  each  new  screen,  two  new  matrices  the  same  size  are  generated.  Where  the 
before  the  memory  was  allocated  for  one  object  of  size  S%,  it  must  now  be  allocated 
for  three  such  objects.  If  N  —  100,  this  means  96  MB  are  now  needed. 

Memory  allocations  are  done  by  segment,  so  when  the  need  exceeds  the  memory 
available,  another  entire  segment  is  allocated.  The  following  operations  use  the  newly 
allocated  memory  until  it  to  is  overfilled,  whereupon,  the  process  of  request  and 
allocation  will  be  repeated  until  there  is  no  more  memory  to  allocate.  However,  the 
request  and  allocation  take  time,  increasing  the  execution  time  of  routines  that  must 
pause  to  do  so.  It  is  this  increase  that  is  reflected  in  the  peaks  in  figure  4.2. 

As  far  as  per-screen  generation  time,  Ts_fs  both  starts  out  higher  and  rises 
more  rapidly  than  both  Ts_z  and  Ts_fft ,  making  the  FS  method  the  least  per-screen 
time  efficient.  Ts_z  is  more  comparable  to  Ts_fft ,  but  the  FFT  method  is  the  most 
per-screen  time  efficient. 

4-1-3  FS  Time  Enhancement.  The  FS  method  can  be  made  to  execute 
more  rapidly  by  decreasing  the  number  of  frequency  components,  controlled  by  the 
multiplicative  frequency  spacing  parameter  Q.  However,  there  a  tradeoff  between 
execution  time  and  number  of  frequencies,  which  becomes  lower  as  Q  increases.  For 
instance,  increasing  Q  from  1.5  to  2.0  retains  76%  of  the  frequencies  while  using  only 
56%  of  the  time:  a  24%  loss  for  a  44%  gain.  But  when  you  increase  Q  from  2.0  to  2.5 
there  is  a  20%  loss  for  a  20%  gain,  and  that  is  about  as  good  as  it  gets  from  there  on 
out. 
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Figure  4.4:  Mean  Time  taken  to  generate  a  100  x  100  pixel 
FS  phase  screen,  with  respect  to  the  Q  value.  Q  is  inversely 
proportional  to  the  number  of  frequency  components. 

The  frequency  grid  remains  polar  and  logarithmic,  so  both  low  and  high  frequen¬ 
cies  are  retained.  They  are  simply  more  sparse  as  Q  increases.  However,  eventually 
the  frequencies  will  become  so  sparse  that  one  loses  the  advantage  gained  by  the  polar 
frequency  grid.  There  are  other  possible  options  for  optimizing  the  FS  method  that 
should  be  explored  as  well.  They  will  be  discussed  in  chapter  5. 

4-1-4  Zernike  Time  Enhancement.  The  Zernike  method  will  also  execute 
faster  if  you  decrease  the  number  of  Zernike  polynomials  used  Nz-  The  difference 
being  that  what  is  lost  with  each  polynomial  is  not  a  specific  frequency  component, 
but  a  particularly  shaped  abberation.  Zernike  polynomials  can  be  selected  to  include 
those  abberations  most  important  to  the  user.  This  may  allow  an  effective  model  to 
be  developed  using  a  very  few  polynomials  that  would  be  extremely  time  efficient. 

Figures  4.5  and  4.6  show  the  initialization  and  per-screen  generation  times, 
respectively,  for  various  numbers  of  Zernike  polynomials.  It  can  be  seen  that  Tj_^ 
increases  quadratically  as  a  function  of  Nz ,  while  Ts_z  increases  linearly.  The  total 
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Figure  4.5:  Time  spent  initializing  100  x  100  pixel  Zernike 

phase  screens  with  respect  to  the  number  of  Zernike  polynomials 
used. 

time  is  second  order  function  of  Nz  as 


Tt~z  —  Ti_z  +  NscrnsTs_z  oc  IVf  +  Nz 


As  the  number  of  phase  screens  Nscrns  generated  increases,  the  quadratic  initialization 
time  factor  becomes  less  dominant. 

An  advantage  that  the  Zernike  method  has  over  the  Fourier  methods  is  that  it 
does  not  require  the  generation  of  new  polynomials  for  changing  atmospheric  condi¬ 
tions.  Such  variations  are  included  in  the  random  weighting  factors  correlated  by  the 
D/r^  proportionate  Zernike  covariance  matrix.  Therefore,  regardless  of  how  many  or 
how  different  the  phase  screens  needed,  a  large  time  investment  in  the  beginning  can 
yield  efficiency  later  on. 
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Figure  4.6:  Time  after  Tj  spent  per-screen  generating  100  x 
100  pixel  Zernike  phase  screens  with  respect  to  the  number  of 
Zernike  polynomials  used.  Does  not  include  T, . 


4-2  Zernike  vs  Fourier  Frequency  Content 

The  Fourier  The  Zernike  polynomials  have  frequency  content  that  tends  to  be¬ 
come  higher  as  the  index  increases,  but  there  is  no  direct  correspondence  to  the 
frequency  content  of  the  FFT  and  FS  sinusoids.  Figure  4.7  shows  a  set  of  Zernike 
polynomials  selected  for  their  increasing  radial  orders  n  and  figure  4.8  shows  their 
Fourier  transforms,  i.e.  their  frequency  content.  Instead  of  each  polynomial  being 
represented  by  a  discrete  spike  in  the  frequency  domain,  as  are  the  sinusoids  from 
the  FFT  and  FS  methods,  they  form  a  continuous  spread,  showing  that  each  one 
is  composed  of  multiple  frequency  components.  However,  there  is  increasing  power 
contained  in  the  higher  frequencies  as  the  Zernike  index  increases. 


4-3  Structure  Function 

The  structure  function  of  an  atmospheric  random  phase  screen,  being  a  mean 
square  difference  function,  should  increase  with  separation  distance,  indicating  that 
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Figure  4.7:  Various  Zernike  polynomials  from  Z2  to  Z74 


i  =  40 


i  =  13 


i  =  67 


Figure  4.8:  The  Fourier  transforms  of  the  zernike  polynomials 
from  Z2  to  Z74  shown  in  figure  4.7,  raised  to  the  1/3  power. 
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Figure  4.9:  Zernike  structure  functions  taken  of  screens  con¬ 
structed  using  Zernike  polynomials  that  form  complete  sets  of 
each  radial  order  n.  From  left  to  right,  and  up  to  down,  Nz  =  2, 

35,  60,  and  84,  respectively.  D/rO  =  1  and  N  =  100. 

correlation  decreases  with  separation  distance  R.  The  slope  of  a  theoretical  structure 
function  on  a  loglog  plot  is  5/3  because  Dn  oc  ( D/r0 )5//3,  where  larger  R  corresponds  to 
smaller  spatial  frequencies  /.  Also,  the  structure  functions  taken  along  the  diagonals 
of  a  phase  map  should  match  those  taken  along  the  axes.  This  shows  the  radial 
symmetry  that  results  from  the  assumption  of  homogeneity  and  isotropy. 

4-3.1  Zernike  Structure  Function.  The  structure  function  of  Zernike  phase 
screens  is  a  not  only  a  function  of  separation  distance  R ,  but  also  of  the  number  and 
radial  order  of  the  Zernike  polynomials  used  to  form  the  phase  screen.  Figures  4.9 
and  4.10  both  show  the  structure  functions  of  four  different  phase  screens  constructed 
for  D/rO  =  1  using  increasing  numbers  of  Zernike  polynomials. 

In  figure  4.9,  all  of  the  screens  were  constructed  using  complete  sets  of  radial 
order  Zernike  polynomials.  The  2  polynomials  used  in  the  upper  left  plot  are  tip 
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and  tilt,  the  complete  set  of  first  radial  order  [n  =  1)  Zernike  polynomials;  the  35 
polynomials  used  in  the  upper  right  plot  form  the  complete  sets  of  Zernike  polynomials 
for  radial  orders  n  —  1  through  n  —  7,  etc.  This  ensures  that  each  phase  screen  will 
be  radially  symmetric,  as  fits  the  assumption  of  homogeneity  and  isotropy. 

The  radial  symmetry  of  a  phase  screen  is  reflected  in  the  similarity  of  the  axial 
structure  function  to  the  the  diagonal  structure  function.  Using  Zernike  polynomials, 
this  symmetry  is  seen  when  n  >  4,  as,  even  given  complete  radial  order  sets,  a  certain 
number  are  required  to  ensure  this  symmetry.  Figure  4.9  gives  an  example:  on  the 
upper  left  only  the  set  of  n  —  1  is  used,  on  the  upper  right  sets  n  =  1  through  n  =  7 
are  used.  One  can  see  that  the  lines  are  almost  entirely  overlapping  in  the  latter  plot. 
The  lower  plots  include  even  greater  numbers  of  radial  order  sets.  In  the  bottom  left 
plot  especially,  it  is  almost  impossible  to  tell  the  axial  structure  function  from  the 
diagonal  structure  function,  indicating  an  almost  exact  radial  symmetry. 

However,  if  the  Zernike  polynomials  are  not  taken  in  radial  order  sets,  radial 
symmetry  can  be  lost,  regardless  of  the  number  of  Zernike  polynomials  used.  Figure 
4.10  illustrates  this  nicely.  In  the  upper  left  is  the  structure  function  of  a  screen  built 
with  only  36  Zernike  polynomials.  On  the  lower  left  is  the  structure  function  of  one 
built  with  60  Zernike  polynomials.  Both  screens  were  built  using  complete  n  sets, 
and  both  have  axial  and  diagonal  structure  functions  that  almost  entirely  overlap. 
On  the  right  are  structure  functions  from  screens  that  each  include  a  partial  n  set. 
The  screens  were  built  with  59  and  61  polynomials,  for  the  upper  and  lower  right 
plots,  respectively,  but  they  are  less  radially  symmetric  than  the  screen  on  the  upper 
left  built  with  only  36  polynomials. 

Increasing  the  number  of  Zernike  polynomials  used,  however,  does  improve  the 
performance  at  higher  frequencies.  One  can  see  in  figure  4.9  that  as  imax  increases, 
high-frequency  correspondence  of  the  Zernike  structure  function  to  the  theoretical 
structure  function.  Hypothetically,  by  including  an  infinite  number  of  Zernike  poly¬ 
nomials,  the  Zernike  phase  screen  structure  function  would  exactly  match  the  theo- 
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Figure  4.10:  Zernike  structure  functions  taken  of  screens  con¬ 
structed  using  Zernike  polynomials  that  do  (on  the  left)  and  do 
not  (on  the  right)  form  complete  sets  of  each  radial  order  n. 
From  left  to  right,  and  up  to  down,  Nz  =  35,  59,  60,  and  61, 
respectively.  D/rO  =  1  and  N  =  100. 
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retical.  It  can  come  close  for  finite  numbers,  but  at  the  cost  of  increased  generation 
time  proportional  to  N §. 

4-3.2  FFT  Structure  Function.  As  can  be  seen  in  the  upper  left  plot  of  figure 
4.11,  the  structure  function  of  an  FFT  phase  screen  corresponds  well  with  theory  for 
small  R  (high  frequencies),  but  falls  of  at  large  R  (low  frequencies).  The  degradation 
of  Dn_FFT  at  large  R  is  due  the  high  correlation  of  opposing  screen  edges.  This 
correlation  is  caused  by  the  lack  of  low  frequency  content  in  the  component  sinusoids 
The  lowest  non-zero  frequency  component  has  a  period  of  exactly  the  screen  width. 

The  more  rapid  degradation  of  the  structure  function  along  the  diagonals  (seen 
best  in  figure  4.13)  is  explained  by  the  larger  (by  a/2)  diagonal  frequency  intervals 
than  axial  frequency  intervals.  The  decrease  in  size  of  the  longest  sinusoidal  period 
increases  the  rapidity  of  noticeable  period  correlation.  This  phenomenon  is  related 
to  aliasing,  and  can  be  compensated  for.  As  is  shown  in  figure  4.11,  by  generating 
larger  than  needed  phase  screens  and  discarding  all  but  the  center  portion  needed  for 
simulation,  the  low  frequency  content  of  an  FFT  phase  screen  can  be  increased. 

If  the  generated  screen  size  is  twice  the  size  of  the  outer  scale,  then  the  lower 
frequencies  will  be  adequately  populated  and  the  structure  function  will  parallel  the¬ 
ory.  This  can  be  related  to  Nyquist  sampling  in  signal  processing.  Given  the  inverse 
relationship  of  frequency  to  distance,  the  longest  separations  of  a  phase  screen  corre¬ 
spond  to  the  lowest  frequencies.  In  order  to  sample  the  lowest  theoretical  frequency, 
one  must  have  a  screen  size  of  twice  the  theoretical  limits  of  correlation. 

For  example,  given  a  necessary  screen  size  of  D  =  1  nr  with  a  sample  size  of 
Ax  =  0.01  m,  and  a  generated  screen  size  of  N  =  D/dx  =  100,  the  lowest  frequency 
included  in  the  screen  will  be 
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Figure  4.11:  FFT  structure  functions  for  the  center  100  x 

100  pixels  of  screens  generated  using  N  =  100,193,448  and 
704  pixels.  Generating  screens  of  excess  size  includes  more  low 
frequency  content. 

If  the  number  of  pixels  (i.e.  samples)  N  is  increased  so  that  N Ax  =  2 L0  where  L0  is 
the  upper  bound  of  the  inertial  subrange,  and  thus  the  outer  limit  of  correlated  IOR, 
then  the  lowest  frequency  represented  becomes 
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As  any  frequencies  lower  than  represent  separations  larger  than  are  theoretically 
correlated  in  their  IOR,  they  should  not  be  present  in  the  phase  screen. 

The  center  100  x  100  pixels  of  a  2Lo  x  2Lo  FFT  phase  screen  can  be  taken  to  form 
the  desired  phase  screen  of  D  =  1  nr,  which  will  have  a  structure  function  exactly 
matching  theory.  Unfortunately,  this  increase  in  quality  requires  a  corresponding 
increase  in  computational  inefficiency  on  the  order  of  NlogN ,  which  presents  as  an 
increase  in  generation  time. 
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Figure  4.12:  FS  structure  functions  taken  of  screens  con¬ 

structed  using  various  Q  values.  There  is  a  graceful  degrada¬ 
tion  of  quality  as  Q  increases,  due  to  the  increasingly  sparse 
frequency  grids. 


4-3.3  FS  Structure  Function.  The  FS  structure  function  parallels  theory 
across  the  spectrum,  though  it  does  degrade  slowly  with  the  number  of  frequencies. 
Figure  4.12  shows  the  structure  functions  of  screens  constructed  using  various  Q 
values,  corresponding  to  numbers  of  frequencies  Nk  from  640  on  the  upper  right 
(Q  =  1.5)  to  120  on  the  lower  left  (Q  =  5.5).  Note  that,  unlike  either  the  FFT  or 
the  Zernke  structure  functions,  the  FFT  structure  function  does  not  fall  off  for  any 
extremity  of  frequency. 

The  variance  from  theory,  presenting  as  ‘waves’  in  the  structure  function,  for 
higher  Q' s  is  due  to  the  sparse  frequency  content.  There  are  simply  not  enough 
sinusoids  to  give  an  accurate  representation  of  random  phases.  Frequencies  that  are 
not  well  represented  show  up  as  peaks  in  the  structure  function  difference  from  theory. 
However,  at  no  frequency  band  is  this  degradation  catastrophic;  all  frequencies  are 
somewhat  represented,  even  using  only  120  sinusoids,  some  just  more  so  than  others. 
This  is  a  benefit  of  the  logarithmic  spacing  of  the  frequency  grid. 
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Figure  4.13:  The  average  structure  functions  for  FFT,  FS  and 
Zernike  phase  screens,  calculated  along  the  axes  (on  the  left)  and 
along  the  diagonals  (on  the  right);  and  the  theoretical  structure 
function  (denoted  Th)  for  comparison.  Parameters:  D/tq  =  1, 

Q  =  2.0,  Nz  =  84,  N  =  100. 

The  FS  axial  and  diagonal  structure  functions  almost  exactly  overlap,  showing 
a  high  degree  of  radial  symmetry  that  is  only  slightly  affected  by  increasing  Q.  This  is 
a  benefit  of  the  polar  frequency  grid.  The  radial  symmetry  and  graceful  degradation 
make  the  FS  structure  function  useful  even  at  high  Q  values.  However,  for  lower  Q 
values  the  FS  structure  function  almost  exactly  matches  theory,  differing  only  by  a 
constant  that  can  be  easily  accounted  for  in  code.  This  makes  it  highly  applicable  to 
all  types  of  simulations. 

4-4  Discussion 

Figure  4.13  shows  the  structure  functions  of  the  three  types  of  phase  screens 
discussed  herein  calculated  along  the  axes  and  the  diagonals,  respectively.  It  also 
shows  the  theoretical  structure  function  for  comparison. 
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The  irregularities  (presenting  as  wavering  near  R  =  1)  of  the  diagonal  of  the 
FFT  and  FS  structure  functions  are  present  because  the  sample  set  for  the  mean 
difference  drops  below  acceptable  statistical  parameters.  In  fact,  the  largest  three 
separations  take  the  mean  square  difference  of  two  matrices  composed  of  only  9,  4, 
and  1  pixels  (for  Rmax  —  3,  2,  and  1,  respectively),  because  it  is  only  the  far  corners 
of  the  screen  that  are  overlapped.  This  does  not  happen  on  the  axial  calculation 
because,  even  at  the  farthest  separation,  there  is  still  a  set  of  N  differences  to  be 
averaged  over. 

4-4- 1  Maximum  Phase  Screen  Size.  It  should  be  noted  that  any  phase  screen 
generated  at  a  size  in  excess  of  L0  should  have  a  structure  function  that  plateaus  at 
R  >  Lq  due  to  the  uncorrelated  random  behavior  of  the  atmospheric  IOR  beyond  that 
point.  I.e.  the  phases  at  separations  greater  than  L0  are  completely  independent,  so 
their  mean  squared  expected  difference  will  reach  a  constant  maximum. 


4-4-2  FFT  vs  FS  Comparison.  It  is  interesting  to  note  that  the  FS  method 
achieves  better  results  than  the  FFT  method  while  using  far  fewer  frequencies.  In 
figure  4.13,  the  FS  method  uses  only  416  frequencies  and  has  a  structure  function 
parallel  to  theory,  while  the  FFT  method  uses  N 2  =  10,  000  frequencies  and  still  falls 
off  at  higher  R ,  failing  to  model  low  frequencies.  However,  more  time  is  required  using 
the  FS  method. 


For  a  1  nr  (N  =  100)  screen  with  a  conservative  upper  boundary  Lq  =  10,  a 
sampling  width  of  Ax  =  0.01  nr,  the  FFT  method  would  have  to  generate  a  NgenxNgen 
pixel  screen  where 


Nnpn  =  —  =  2,  000 


9en 


Ax 


and  Ngen  =  4  x  106  is  the  total  number  of  pixels  per  screen.  This  would  require  a 


number  of  operations  on  the  order  of  Ngen  log (Ngen)  m  6,  600  and  a  corresponding 
generation  time.  Extrapolating  from  figures  4.1  and  4.3,  this  leads  to  a  total  genera- 
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tion  time  for  two  phase  screens  on  the  order  of  3  to  5  seconds,  the  majority  of  which 
is  the  per-screen  generation  time  Ts. 

On  the  other  hand,  the  FS  method,  generating  a  screen  of  similar  quality  for 
the  same  specifications,  would  require  only  Nk  =  416  frequencies  (Q  =  2.0).  The  FS 
screen  generated  would  be  exactly  the  size  needed  so  the  number  of  operations  would 
be  on  the  order  of  N2Nx  =  4.16  x  106  due  to  the  lack  of  optimization  within  the  code. 
From  figures  4.1  and  4.2,  it  takes  6  to  7  seconds  to  generate  two  equivalent  100  x  100 
phase  screens,  most  of  which  constitutes  the  initialization  time  T*. 

Because  most  of  the  time  for  the  FS  screen  is  spent  building  the  kernel,  sub¬ 
sequent  phase  screens  can  be  generated  at  the  rate  of  0.5  second  per  pair.  But 
because  most  of  the  FFT  generation  time  is  spent  generating  the  screen  itself,  addi¬ 
tional  screens  will  take  3  seconds  per  pair.  Therefore,  in  a  situation  where  quality  is 
important  and  more  than  four  phase  screens  are  required,  it  is  more  time  efficient  to 
use  the  FS  method  than  the  FFT  method. 
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V.  Conclusions  and  Recommendations 


These  conclusions  and  recommendations  stem  directly  from  the  results  covered 
in  chapter  4.  They  include  suggested  applications  of  the  various  methods  as 
well  as  possible  optimizations  that  might  be  made. 

5. 1  Zernike  Method 

The  Zernike  is  highly  effective  at  capturing  low  frequency  abberations,  such  as 
tip  and  tilt.  It  is  also  radially  symmetric  for  complete  radial  order  sets  of  Zernike 
polynomials  beyond  n  >  4.  This  makes  it  highly  effective  for  simulating  turbulence 
when  low  frequency  accuracy  is  required,  such  as  for  use  with  tracking  and  target¬ 
ing  algorithms.  The  main  fault  with  the  Zernike  method  is  that  it  falls  off  at  high 
frequencies  for  limited  numbers  of  polynomials;  however,  including  the  high  numbers 
(on  the  order  of  500  to  1000)  of  polynomials  needed  to  accurately  model  the  high 
frequencies  can  become  cumbersome  as  far  as  generation  time  is  concerned. 

5.2  FFT  Method 

The  FFT  method  is  by  far  the  fastest  for  generating  phase  screens,  but  it  di¬ 
verges  from  theory  at  low  frequencies  for  screen  sizes  of  less  than  2 L0,  and  becomes 
cumbersome  at  sizes  near  2Lq.  This  makes  it  good  for  simulating  high  frequency  tur¬ 
bulence  for  adaptive  optics  algorithms,  where  low  frequency  abberations  are  assumed 
to  be  already  removed. 

5.3  FS  Method 

The  FS  method  takes  the  longest  initialization  time,  but  it  has  the  best  results 
across  the  spectrum.  It  does  not  fail  catastrophically  at  any  frequency  range,  instead 
it  degrades  gracefully  with  few  frequencies  modelled.  FS  phase  screens  are  applicable 
to  any  scenario,  but  are  computationally  cumbersome  as  implemented  here.  The  FS 
method  is  therefore  well  worth  optimizing.  One  obvious  area  for  improvement  is  to 
find  a  way  other  than  a  triply-nested  for  loop  to  construct  the  kernel. 
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Given  an  optimized  implementation,  the  FS  method  would  be  universally  ap¬ 
plicable,  generating  quality  results  across  the  spectrum.  However,  even  as  is,  the  FS 
method  has  a  comparable  implementation  time  with  either  the  Zernike  or  the  FFT 
methods  for  the  same  quality  of  result. 
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Appendix  A.  MATLAB  Code 

Note  that  filenames  have  had  underscores  removed  and  capitals  inserted  for  LaTex 
display  calls.  The  hies  must  be  renamed  to  match  the  code  before  implementation. 


Listing  A.l:  Top  level  code  to  generate  all  three  types  of  phase  screens, 
(appendixl  /  GeneratePhzScrnsDn.m) 


°/0  GeneratePhzScrnsDn.m 
7„  Rebecca  Eckert 
1  15  Feb  06 


°/.  °/.  y.  y.  y.  y.  y  y.  •/.  y.  y  y  y.  %  •/.  y.  y.  y.  y.  y y  y.  y.  %  %  y.  y.  y  y.  y.  y  y  y.  y.  %  %  y.  y  y  y y y  y.  y.  y.  y.  y.  y.  y.  y.  y.  y.  y.  y.  y.  %  y.  y.  y.  y  y.  y.  y.  y.  y.  y.  %  y.  y.  y y y y y y 


°/0  GENERATEPHZSCRNSDN.M  Generates  Zernike  ,  FFT  phase  screens  and 
calculates 

•/.  their  structure  functions  given  a  particualr  D/rO.  It  also  ... 
records  the 

•/.  time  taken  to  generate  each  screen, 
y.  -  Iterate  D/rO  =  [1  3  5  10] 

1  -  Iterate  30  phase  screens/Sturcture  Functions  per  D/rO 


y.  y.  y.  y.  %  %  y.  y.  y  y.  y.  y.  y.  y.  y.  y  %  %  y.  y  y  y  %  y.  y.  y.  y.  y.  %  y.  y.  y  y.  y  y.  y  y.  y.  y  y.  y.  y.  y.  y.  y  y  y  y.  y.  y.  y.  y.  y.  %  %  %  y  y  y. 


y.  y.  y.  y.  y.  %  %  y  y y y y y y 1  y. 


clear ; clc ; matlabpath (pathdef ) 


y  inputs  y.  y.  y.  y.  y.  y  y  y„  y„  y0  y„  y0  y.  y.  y0  yt  y0  yt  y„  y„  y0  y„  y.  y„  y.  y.  y.  y.  y.  y„  y„  y0  y„  y.  %  y0  y0  y0  y,  yt  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y 


y  Screen  Descriptors 
N  =  100; 

NumScreens  =  60; 
generated 

NumZernikeModes  =  84; 
used  in 

•/.  Atmospheric  Descriptors 
D  =  1; 

rO  =  [D]  0/o  D/3  D/5  D/10]  ;% 
AtmSpectrum  =  ’VonKarmen' 
’  ,  ’ ModAtm  ’ 


“/,  Large  inertial  subrange 
lo  =  0.01; 

Lo  =  10; 

0/.  Small  Inertial  subrange 
0/.  lo  =  0.01; 

0/.  Lo  =  10; 


y  The  size  in  pixels  of  phase  screens 
°/o  The  number  of  phase  screens  to  be  ... 

0/o  The  number  of  polynomials /mode  to  be... 

“/o  generating  the  Zernike  Phase  Screen 

“/,  optic  Diameter 
Fried  radius 

“/„  1  Kolomogorov  ’  ,  ’  Tatarski  ’  ,  ’  VonKarmen  .  .  . 

y  =>  Used  in  making  the  FFT  phase  ... 
screen 

0/o  Inertial  Subrange  Inner  Scale  (m) 

“/o  Inertial  Subrange  Outer  Scale  (m) 


lambda  =  0.5e-6;  Wavelength  (m) 

WaveNumber  =  2*pi/lambda;  °/n 


”/o  Reference  Coordinates  -  Spatial  Domain 
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R  =  linspace (0 , D) ; 
dx  =  D/N; 


7. 7 7 7. 7. 7 7 7 7 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7 7 7 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7 


for  kk  =  l:length(rO) 

D_r0  =  D/rO (kk) ; 

PhzFileName  =  [ ’ PhzScrns_DrO_  ’  num2str (D_rO) ] 


7 


Ideal  structure  function... 


7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7 


Dphi_Ideal  =  6 . 88* (R/rO (kk) )  .  ~  (5/3)  ; 


7  FFT  PHASE  SCREEN... 


7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7 


tic 

PlizScrnFFT  =  zeros  (N  ,  N  ,  NumScreens  )  ; 
ind  =  1 ; 

for  ii  =  1 : NumScreens /2 

[PI , P2 ]  =  MakeFFT_PhzScrn(D,rO(kk)  ,Lo,lo,N)  ; 
PhzScrn_FFT(:,:,ind)  =  PI; 

PhzScrn_FFT(: , : , ind+1)  =  P2 ; 
ind  =  ind  +  2 ; 

end 

toe 


7  ZERNIKE  PHASE  SCREEN... 

7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7 

7  Note  that  the  My_get _zernike_mode . m  function  only  contains  the  ... 
first  85 

7  Zernike  polynomial  modes,  starting  with  Piston  [1].  HOWEVER,  ... 
the 

7  Zernike_Cov  matrix  contains  2999x2999  covariance  terms,  starting... 
with  tip 

7  and  tilt  [1,2]  and  ignoring  piston.  Therefore,  only  84  ... 
polynomials/modes 

7  are  currently  available,  limiting  the  higher  frequency  ... 
capabilities . 

7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7 


7  Zernike  Coevariance  Matrix 
load  z_cov . mat 

zernike  =  My_get_zernike_mode (NumZernikeModes  +  1,N); 
zernike  =  zernike (2: NumZernikeModes  +  1  ,  :  ,  :  )  ; 
load  zernike_index 
ZernikeCov  =  z_cov; 

tic 

PhzScrn_Zern  =  zeros (N , N , NumScreens ) ; 
for  ii  =  1 : NumScreens 

PhzScrn_Zern(:,:,ii)  =  ... 

MakeZernikePhzScrn(D_rO , NumZernikeModes  , ZernikeCov ,N,  . .  . 
zernike )  ; 

end 
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toe 


80 


85 


90 


7.  Housekeepint  to  optimize  FS  routine 

save  (  PhzFi  leName  ,  1  PhzScrn  *’,  ’  dx  ’  ,  ’  Dphi  _  Ideal  R  ’  )  °/  Keep  phase 

s  erns 

save  tempi  D  rO  N  kk  lo  Lo  D_rO  NumScreens  PhzFileName 

clear  PhzScrn* 

save  temp2 

clear  all 

load  tempi 

/  RADIAL  FOURIER  SERIES  PHASE  SCREENS... 


7.  / 1 1 1 1 7. 7. 7. 7. 7, 7. 7. 7. 7. 1 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 


7,  Generates  a  screen  similar  to  the  FFT  screen,  but  with 
logr ithmi c 

7.  frequency  components  generated  on  a  polar  grid. 


7. 7. 7. 7. 7. 7. 7. 1 1 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 1 1 1 1 1 1 1 1 1 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 1 1 1 7. 7. 1 1 1 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 


tic 

x_p  =  linspace (0  ,  . 5 , N)  ; 

95  y_p  =  x_p ; 

Q  =  2.0 

kern  =  My_make_polar_kern (x_p , y_p , rO (kk)  , lo , Lo  ,  Q)  ; 

toe 

tic 

100  ind  =  1  ; 

for  i i = 1 : NumScreens /2 

[phil,phi2]  =  make_polar_screen (kern) ; 
PhzScrn_FS(:,:,ind)  =  phi 1 -mean (phil  (:))  ; 
PhzScrn_FS(:,:,ind  +  1)  =  phi2-mean (phi 2 ( : ) ) ; 
105  ind  =  ind  +  2 ; 

end 
toe 


save (PhzFileName , ’-append1 , ’PhzScrn_FS’ , ’ x_p ’ , ’  Q  ’  ) 
110  clear  PhzScrn_FS 


7. 7. 7. 7. 7. 1 1 1 1 1 1 1 1 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 1 1 1 1 1 1 1 1 1 1 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7.  %  7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 7. 


load  temp2 
end 


Listing  A. 2:  Generate  a  Zernike  phase  screen. 

(appendixl/MakeZernikePhzScrn.m) 

function  [  PhzScrn  ]  =  MakeZernikePhzScrn (D_r0 , numModes  , z_cov , . .  . 

numPix , zernike) 

7,  MAKEZERNIKEPHZSCRN  Creates  a  random  phase  screen  using  the  ... 
Zernike 
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7o  technique.  The  phase  screen  is  generated  in  cartesian  ... 

coordinates  on  a 
7o  rectangular  grid. 

7o  NOTE  that  the  Zernike  will  be  a  circle  inscribed  on  the  ... 

rectangular 
7.  frequency  grid. 

7.  NOTE  that  the  My_get  _zernike_mode  .  m  function  only  contains  the 
first  85 

7.  Zernike  polynomial  modes,  starting  with  Piston  [1]. 

7.  Generate  Zernike  Phase  screen 
if  ( isempty (numPix )  ==  1) 

error (’ MakePhzScrn  of  type  =  ‘ZERNIKE1  requires  input  SIZE  in 
pixels’) 

elseif  ( i sempt y ( D_r 0 )  ==  1) 

error ( 1 MakePhzScrn  of  type  =  ‘ZERNIKE1  requires  input  ... 
atmospheric  D/rO’) 
else  npix  =  numPix; 
end 

if  (mod(npix,2)  ==  1) 

[x  y]  =  meshgrid (-floor (npix/2)  : floor (npix/2) )  ; 

else  [x  y]  =  meshgrid (npix/2 : npix/2-1) ; 

end 

r  =  sqrt(x.~2  +  y.~2); 

7.  Zernike  Coevariance  Matrix 
7.  load  z_cov.mat 

7.  zernike  =  My_get_zernike_mode  (numModes  +  l,npix); 

7.  zernike  =  zernike  (2  :  numModes  +  1  ,  :  ,  :  )  ; 

a  =  chol(z_cov)  ’* randn ( length ( z_cov )  ,1)  ; 
a_ps  =  a ( 1 : numModes )* D_rO ~ ( 5/3 ) ; 

Aps  =  ( r epmat ( a_ps , [ 1  npix  npix])); 

PhzScrn  =  squeeze ( sum ( Aps .* zernike  ,  1 ))  ; 


Listing  A. 3;  Generate  FFT  phase  screen. (Using  Von  Karmen  spectral  model.) 

(appendixl/MakeFFTPhzScrn.m) 


function  [  PhzScrnl  PhzScrn2  ]  =  MakeFFT_PhzScrn (D , rO , Lo , lo , N) 

7.  MAKEFFT_PHZSCRN  Creates  a  random  phase  screen  using  the  FFT  ... 
technique 

%  and  the  Von  Karmen  spectral  model 

7. 7. 7. 7. 7. 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 


DrO  =D/rO ; 

[p  q]  =  meshgr id ( ( -N /2 : N/2 - 1 ) ) ; 
rl  =  randn ( N )  +  j*randn(N); 
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10  G  =  (p.~2  +  q.~2  +  (D/Lo) ~2)  . ~ ( -11/12)  . *exp (-0 . 563*(lo/D) “2* (p . “2+q. .  . 

.-2)); 

G (N/2+1 , N/2+1) =0; 

F  =  f f tshif t ( fft2 ( f ft  shift ( G  .  *rl ))) ; 

F  =  F*0 . 1517*  DrO ~  (5/6)  ; 

15  PlizScrnl  =  real(F); 

Ph.zScrn.2  =  imag(F); 


7,  isc  (PhzScrnl  )  ;  colorbar 

20 


%  OLD  CODE 

°/«  if  (  i  sempty  (  PHI  )  ==  1) 

25  7.  error  (  1  MakePhzScrn  of  type  =  ‘  FFT  ‘  ... 

7.  requires  input  PHI  =  atmospheric  spectrum  (eg,  Von  Karmen)  ’) 

7.  elseif  (  isempty  (numPix)  ==  1)  I  (numPix  <  length  (PHI  )  ) 

/  npix  =  length(PHI); 

7.  else 

30  npix  =  length(PHI); 

7.  end 

7.  S  =  sqrt(PHI); 

7.  Generate  Gaussian  random  phase  screen  using  FFT  technique 
g  =  randn(npix); 

35  PhzFFT  =  g  .  *PHI ; 

G  =  f ft  shift ( if ft  2 ( f ft  shift ( PhzFFT )))  ; 

PhzScrnl  =  real(G); 

PhzScrn2  =  imag(G); 

Listing  A. 4:  Generate  FS  phase  screen. (appendixl/MakePolarScreen.m) 

/.  =======  =  =  =  =  ======  =  =  =  =  =  =  =====  =  =  =  =  ======  =  =  =  =  =======  =  =  =  =  =======  =  =  =  = 


7.  Creates  a  pair  of  phase  screens  with  subharmonics  sampled  over  a... 

pseudo  - 
7.  polar  gr  id  . 

/.  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  = 


5  function  [phil,phi2]  =  make_polar_screen (kern) 

[Ny,Nx,Ns]  =  size(kern); 

comp  =  randn  ( 1  , 1  ,  Ns  )  +  j  * r andn  (  1  ,  1  ,  Ns  )  ; 
phi  =  zeros  (Ny  ,  Nx  )  ; 

10  for  il  =  1 : Nx 

for  i2  =  1  :  Ny 

phi(i2,il)  =  phi(i2,il)  +  sum (comp.*kern(i2,il  ,  :)  ,3)  ; 

end 

end 

15  7.  phi_s  =  sum  (kern  .*shiftdim(  repmat  (comp,[l,Ny,  Nx]),l),3); 
phil  =  real(phi); 
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ph.il  =  phil -mean  (mean  (phil  ))  ; 
phi2  =  imag(phi); 

20  phi2  =  phi2 -mean (mean (phi2 )) ; 

return 


Listing  A. 5:  Generate  FS  kernel.  (appendixl/MyMakePolarKern.m) 

function  kern  =  make.polar _ker (x , y , rO , 10 , LO , Q ) 

7.  7.  Q  =  2.0;  %  Factor  to  set  number  of  points 

next_inc  =  10;  70used  to  display  progress 
5  prog_step  =  10;  70  used  to  display  progress 

max.iter  =  1000;  70  num  of  iterations  for  midpoint  selection 

tol  =  le-12;  7,  tolerance  for  midpoint  selection 

10 

dx  =  x(2)  -  x(l); 

dy  =  y  (2)  -  y  (1)  ; 

Nx  =  length (x) ; 

Ny  =  length (y)  ; 

15 

kappaO  =  1 / L0 ; 

7.  Set  max  freq  also  based  on  sampling 
kappaml  =  5.92/10; 
kappam2  =  2*pi/(2*dx); 

20  kappam  =  min (kappaml , kappam2 ) ; 

k_max  =  kappam; 

k_min  =  2*pi/5/L0;  7«  half  the  lowest  frequency  demanded  by  the  . 
outer  scale 

25  disp (’ Forming  polar  grid  kernel’); 

log_k_min  =  loglO (k_min)  ; 
log_k_max  =  loglO (k_max)  ; 
delta_log  =  log_k_max  -  log_k_min; 

30  num_freq_pts  =  ce il ( delt a_log / logl 0 ( Q ) )  ; 

log_k_pts  =  linspace ( log_k_min , log_k_max , num_f req_pts ) ; 
k_pts  =  10  .  ~  log_k_pts  ;  7«  These  are  the  radial  frequency  values 

used 


35  k_low  =  k_pts(l:end  -  1); 

k_low2  =  (k_pts(l:end  -  1)).~2; 
k_high  =  k_pts (2 : end) ; 
k_high2  =  ( k_pt s (2 : end ) )  .  ~ 2 ; 
delta_k  =  diff(k_pts); 

40 

°i  need  to  find  the  2D  power  in  each  annular  segment  from  PSD 
d_phi2  =  zeros ( 1 , length ( delta_k )) ; 
for  il  =  1 : length ( delta_k ) 
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d_phi2(il)  =  int _vK2D2 ( rO , 10 , L0 , k_pt s ( i 1 ) , k_pt s ( i 1  +  1)); 


end 


7. 7. 7. 7. 7. 1 1 1 1 1  loop  to  pick  point  to  use  for  sinusoids 
7.  define  differential  area 
k_area  =  pi*(k_high2  -  k_low2); 

7.  7.  provide  a  guess  at  the  desired  midpoint  and  then  iterate  to  .. 
choose  best 

7.  7.  midpoint  as  defined  by  place  where  power  in  annular  segment  is 
equal  to 

7.  7.  product  of  area  and  PSD  at  that  midpoint 
k_mid  =  0.5*delta_k  +  k_low; 

7.  cnt  =  0 ; 

7.  max.iter  =  1000; 

7.  tol  =  le-12; 

7.  for  index  =  1  :  length  (k_mid) 


7,  seed  binary  searc  with  endpoints  and  midpoint 
k_lower  =  k_pt s ( index ) ; 
k_upperu  =  k_pts (index  +  1)  ; 

B  =  k_mid ( index ) ; 

phi_mid  =  vK_phi ( B ~2 , rO , 10 , L0 )  ; 
ip  =  phi_mid*k_area ( index) ; 

while  (cnt  <=  max.iter)  &&  (abs(ip  -  d_phi2 ( index ) )  >  tol) 

if  ip  <  d_phi2 ( index ) 
k_upper  =  B; 

B  =  k_lower  +  ,5*(B  -  k.lower)  ; 

else 

k_lower  =  B; 

B  =  k.lower  +  ,5*(k_upper  -  B); 

end 

phi_mid  =  vK_phi (B "2 , rO , 10 , L0 ) ; 
ip  =  phi_mid*k_area ( index) ; 
cnt  =  cnt  +  1  ; 
if  cnt  ==  max.iter 

di sp ( 1  Warning  maximum  iterations  reached  in  kappa 


calc  ’  )  ; 
))]) 


di sp  (  [  ’ 


tol:  ’  num2str ( abs ( ip  -  d_phi2 ( index ) 


end 

end 

k_mid ( index) 
cnt  =  0; 


7.  end 


7.  this  can  be  changed  later  to  include  a  progression  of  theta 
slices 

num.theta  =  32;  %  number  of  slices  in  theta 

theta  =  linspace (pi/num_theta  ,  2*pi  -  pi/num_theta , num.theta)  ; 
theta  =  repmat (theta  ,  [1 , length (k_mid) ])  ; 

k  =  []  ; 
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d_area  =  []  ; 

for  il  =  1 : length (k_mid) 

k  =  [k  repmat(k_mid(il),[l, num_ theta] ) ] ; 

d_area  =  [d_area  r epmat ( k_area ( i 1 ) /num_thet a , [1 , num_thet a] ) ] ; 

end 

kx  =  k  .*  cos  (theta)  ;  "/,  cartesian  points  for  given  theta  value 
ky  =  k .* sin (theta) ; 

sqrt_phi  =  sqrt ( vK_phi (k . ~2 , rO , 10 , LO ) ) ; 
sqrt.area  =  sqrt (2*pi*d_area) ; 

tot_iter  =  Nx*Ny* length (k) ; 
kern  =  zeros (Nx , Ny , length (k) ) ; 
for  il  =  1 : Nx 

for  i2  =  1  :  Ny 

for  i3  =  l:length(k) 

kern ( i2  ,  i 1  ,  i3 )  =  sqrt _phi ( i3 )* sqrt .area ( i3 )*..  . 
exp(j*(x(il)*kx(i3)  +  y ( i2 ) *ky ( i3 ) ) ) ; 

end 

if  i 1  * i2 * i3 *  100/ t ot _ it er  >  next.inc 
disp (num2str (next_inc)) ; drawnow 
next_inc  =  next.inc  +  prog_step; 

end 

end 

end 

disp ( ’ done  !  ’  )  ; 

/  test  code 

7.  output  the  resulting  set  of  kappa  values  for  analysis 
ktst2  =  logspace (loglO (k_pts (1) )  , loglO (k_pts (end) )  ,  1000)  ; 
ptst  =  vK_phi (ktst2 . ~2 , rO , 10 , L0)  ; 

7o  7.  figure  (2)  ;  elf  ; 

7.  7.  subplot  (  1 , 2 , 1 )  ; 

7.  7.  loglog  (ktst2  ,  ptst  ,  1  g- .  ’  )  ; 

7.  7.  hold  on; 

7.  7.  loglog  (k_pts  ,  vK_phi  (k_pts  .  ~2  ,  rO  ,  10  ,  L0)  ,  ’  k+  ’); 

7.  /  loglog  (k_mid  ,  vK_phi  (k_mid  .  "2  ,  rO  ,  10  ,  L0)  ,  1  ro  ’); 

7.  7.  axis  tight; 

7.  /  subplot  (  1  ,  2 , 2)  ; 

7.  7.  semilogy(ktst2,ptst,,g-.’); 

7.  /  hold  on; 

7.  7.  semilogy(k_pts,vK_phi(k_pts.~2,r0,10,L0),,k+  ’ )  ; 

7.  7.  semilogy(k_mid,vK_phi(k_mid.~2,r0,10,L0),,ro  '); 

7.  7.  axis  tight; 

7.  7.  hold  off  ; 

kmtstsq  =  sqrt(k_mid); 

ptstsq  =  sqrt ( vK_phi (k_mid . ~2 , rO , 10 , L0) ) ; 

7.  end  test  code 
return 
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Listing  A. 6:  Calculate  Von  Karmen  spectral  model. (appendixl/vKphi.m) 

7,  calculate  the  von  Karman  phase  spectrum  for  a  given  kappa,  rO  ,  . 

10 ,  and  LO 

function  PHI  =  vK_phi ( kappa2 , rO , 10 , LO ) 

7,  kappa_max  =  5.92/10; 

°l  kappa_max2  =  kappa_max~2; 

5  7.  kappa_min  =  1/LO; 

7.  kappa_min2  =  kappa_min~2; 

7.  PHI  =  0.4916693*r0~(-5/3)*exp(-kappa2./kappa_max2).*... 

7.  (kappa2  +  kappa_min2  )  .  “  (  - 1 1  /6 )  ; 

10  PHI  =  0.4916693*r0~(-5/3)*exp(-kappa2  .  / ( 5 . 92/10)  ,~2)  .*..  . 

( kappa2  +  (1/LO) .“2) .“(—11/6) ; 
return 
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